Multiple tagging of individual long dna fragments

ABSTRACT

This disclosure provides methods and compositions for tagging long fragments of a target nucleic acid for sequencing and analyzing the resulting sequence information in order to reduce errors and perform haplotype phasing, for example.

RELATED APPLICATION

This disclosure is a continuation of U.S. patent application Ser. No.15/136,780, filed Apr. 22, 2016, which is a continuation of U.S. patentapplication Ser. No. 14/205,145, filed Mar. 11, 2014, now U.S. Pat. No.9,328,382, which claims the priority benefit of U.S. provisional patentapplication 61/801,052, filed Mar. 15, 2013. The priority applicationsare hereby incorporated herein by reference in their entirety for allpurposes.

REFERENCE TO A “SEQUENCE LISTING,” A TABLE, OR A COMPUTER PROGRAMLISTING APPENDIX SUBMITTED AS AN ASCII TEXT FILE

The Sequence Listing written in file 92171-902410_ST25.TXT, created May20, 2014, 770 bytes, machine format IBM-PC, MS-Windows operating system,is hereby incorporated by reference.

BACKGROUND OF THE INVENTION

There is a need for improved methods for determining the parentalcontribution to the genomes of higher organisms, i.e., haplotype phasingof genomes. Methods for haplotype phasing, including computationalmethods and experimental phasing, are reviewed in Browning and Browning,Nature Reviews Genetics 12:703-7014, 2011.

Most mammals, including humans, are diploid, with half of the homologouschromosomes being derived from each parent. Many plants have genomesthat are polyploid. For example, wheat (Triticum spp.) have a ploidyranging from diploid (Einkorn wheat) to quadriploid (emmer and durumwheat) to hexaploid (spelt wheat and common wheat [T. aestivum]).

The context in which variations occur on each individual chromosome canhave profound effects on the expression and regulation of genes andother transcribed regions of the genome. Further, determining if twopotentially detrimental mutations occur within one or both alleles of agene is of paramount clinical importance. For plant species, knowledgeof the parental genetic contribution is important for breeding progenywith desirable traits.

Current methods for whole-genome sequencing lack the ability toseparately assemble parental chromosomes in a cost-effective way anddescribe the context (haplotypes) in which variations co-occur.Simulation experiments show that chromosome-level haplotyping requiresallele linkage information across a range of at least 70-100 kb. Thiscannot be achieved with existing technologies that use amplified DNA,which are be limited to reads less than 1000 bases due to difficultiesin uniform amplification of long DNA molecules and loss of linkageinformation in sequencing. Mate-pair technologies can provide anequivalent to the extended read length but are limited to less than 10kb due to inefficiencies in making such DNA libraries (due to thedifficulty of circularizing DNA longer than a few kb in length). Thisapproach also needs extreme read coverage to link all heterozygotes.

Single molecule sequencing of greater than 100 kb DNA fragments would beuseful for haplotyping if processing such long molecules were feasible,if the accuracy of single molecule sequencing were high, anddetection/instrument costs were low. This is very difficult to achieveon short molecules with high yield, let alone on 100 kb fragments.

Most recent human genome sequencing has been performed on shortread-length (<200 bp), highly parallelized systems starting withhundreds of nanograms of DNA. These technologies are excellent atgenerating large volumes of data quickly and economically.Unfortunately, short reads, often paired with small mate-gap sizes (500bp-10 kb), eliminate most SNP phase information beyond a few kilobases(McKernan et al., Genome Res. 19:1527, 2009). Furthermore, it is verydifficult to maintain long DNA fragments in multiple processing stepswithout fragmenting as a result of shearing.

At the present time three personal genomes, those of J. Craig Venter(Levy et al., PLoS Biol. 5:e254, 2007), a Gujarati Indian (HapMap sampleNA20847; Kitzman et al., Nat. Biotechnol. 29:59, 2011), and twoEuropeans (Max Planck One [MP1]; Suk et al., Genome Res., 2011;http://genome.cshlp.org/content/early/2011/09/02/gr.125047.111.full.pdf;and HapMap Sample NA 12878; Duitama et al., Nucl. Acids Res.40:2041-2053, 2012) have been sequenced and assembled as diploid. Allhave involved cloning long DNA fragments into constructs in a processsimilar to the bacterial artificial chromosome (BAC) sequencing usedduring construction of the human reference genome (Venter et al.,Science 291:1304, 2001; Lander et al., Nature 409:860, 2001). Whilethese processes generate long phased contigs (N50s of 350 kb [Levy etal., PLoS Biol. 5:e254, 2007], 386 kb [Kitzman et al., Nat. Biotechnol.29:59-63, 2011] and 1 Mb [Suk et al., Genome Res. 21:1672-1685, 2011])they require a large amount of initial DNA, extensive libraryprocessing, and are too expensive to use in a routine clinicalenvironment.

Additionally, whole chromosome haplotyping has been demonstrated throughdirect isolation of metaphase chromosomes (Zhang et al., Nat. Genet.38:382-387, 2006; Ma et al., Nat. Methods 7:299-301, 2010; Fan et al.,Nat. Biotechnol. 29:51-57, 2011; Yang et al., Proc. Natl. Acad. Sci. USA108:12-17, 2011). These methods are useful for long-range haplotypingbut have yet to be used for whole-genome sequencing; they requirepreparation and isolation of whole metaphase chromosomes, which can bechallenging for some clinical samples.

There is also a need for improved methods for obtaining sequenceinformation from mixtures of organisms such as in metagenomics (e.g.,gut bacteria or other microbiomes). There is also a need for improvedmethods for genome sequencing and assembly, including de novo assemblywith no or minimal use of a reference sequence), or assembly of genomesthat include various types of repeat sequences, including resolution ofpseudogenes, copy number variations and structural variations,especially in cancer geneomes.

We have described long fragment read (LFR) methods that provide enablean accurate assembly of separate sequences of parental chromosomes(i.e., complete haplotyping) in diploid genomes at significantly reducedexperimental and computational costs and without cloning into vectorsand cell-based replication. LFR is based on the physical separation oflong fragments of genomic DNA (or other nucleic acids) across manydifferent aliquots such that there is a low probability of any givenregion of the genome of both the maternal and paternal component beingrepresented in the same aliquot. By placing a unique identifier in eachaliquot and analyzing many aliquots in the aggregate, DNA sequence datacan be assembled into a diploid genome, e.g., the sequence of eachparental chromosome can be determined. LFR does not require cloningfragments of a complex nucleic acid into a vector, as in haplotypingapproaches using large-fragment (e.g., BAC) libraries. Nor does LFRrequire direct isolation of individual chromosomes of an organism. Inaddition, LFR can be performed on an individual organism and does notrequire a population of the organism in order to accomplish haplotypephasing. LFR methods have been described in U.S. patent application Ser.Nos. 12/329,365 and 13/447,087, U.S. Pat. Publications US 2011-0033854and 2009-0176234, and U.S. Pat. Nos. 7,901,890, 7,897,344, 7,906,285,7,901,891, and 7,709,197, all of which are hereby incorporated byreference in their entirety.

BRIEF SUMMARY

The present invention provides methods and compositions for MultipleTagging of Individual Long DNA Fragments (referred to herein by theabbreviation Multiple Tagging, or MD and sequencing and for analyzingthe resulting sequence information in order to reduce errors and performhaplotype phasing, among other things.

According to one aspect of the invention, methods are provided forsequencing a target nucleic acid comprising: (a) combining in a singlereaction vessel (i) a plurality of long fragments of the target nucleicacid, and (ii) a population of polynucleotides, wherein eachpolynucleotide comprises a tag and a majority of the polynucleotidescomprise a different tag; (b) introducing into a majority of the longfragments tag-containing sequences from said population ofpolynucleotides to produced tagged long fragments, wherein each of thetagged long fragments comprises a plurality of the tag-containingsequences at a selected average spacing, and each tag-containingsequence comprises a tag; and (c) producing a plurality of subfragmentsfrom each tagged long fragment, wherein each subfragment comprises oneor more tags. Such methods are suitable for preparing a target nucleicacid for nucleic acid sequencing. According to one embodiment, suchmethods further comprise sequencing the subfragments to produce aplurality of sequence reads; assigning a majority of the sequence readto corresponding long fragments; and assembling the sequence reads toproduce an assembled sequence of the target nucleic acid.

According to one embodiment, producing the subfragments in such methodscomprises performing an amplification reaction to produce a plurality ofamplicons from each long fragment. According to another embodiment, eachamplicon comprises a tag from each of the adjacent introduced sequencesand a region of the long fragment between the adjacent introducedsequences.

According to another embodiment, such methods comprise combining thelong fragments with an excess of the population of tag-containingsequences.

According to another embodiment, such methods comprise combining thelong fragments with the tag-containing solution under conditions thatare suitable for introduction of a single tag-containing sequence into amajority of the long fragments.

According to another embodiment, such methods comprise combining thelong fragments with the tag-containing solution under conditions thatare suitable for introduction of a different tag-containing sequencesinto a majority of the long fragments.

According to another embodiment, in such methods the population oftag-containing sequences comprises a population of beads, wherein eachbead comprises multiple copies of a single tag-containing sequence.

According to another embodiment, in such methods the population oftag-containing sequences comprises a population of concatemers, eachconcatemer comprising multiple copies of a single tag-containingsequence.

According to another embodiment, in such methods the tag-containingsequences comprise transposon ends, the method comprising combining thelong fragments and the tag-containing sequences under conditions thatare suitable for transposition of the tag-containing sequences into eachof the long fragments.

According to another embodiment, in such methods the tag-containingsequences comprise a hairpin sequence.

According to another embodiment, in such methods the target nucleic acidis a complex nucleic acid.

According to another embodiment, in such methods the target nucleic acidis a genome of an organism.

According to another embodiment, such methods comprise determining ahaplotype of the genome.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a method for tagging and fragmenting of long fragments of atarget nucleic acid with transposon-mediated barcodes.

FIGS. 2A and 2B show a method for tagging and fragmenting of longfragments of a target nucleic acid with hairpin-mediated barcodes.

FIG. 3 shows a method for transposon-mediated tagging and fragmenting oflong fragments of a target nucleic acid.

FIG. 4A shows a method for tagging long fragments of a target nucleicacid using a tagged adaptor.

FIG. 4B shows an alternative method for tagging long fragments of atarget nucleic acid using a tagged adaptor.

FIGS. 4C-D show a second alternative method for tagging long fragmentsof a target nucleic acid using a tagged adaptor.

FIGS. 4E-F show methods for creating a series of tagged subfragments ofwith shorter and shorter regions of the long DNA fragments.

FIGS. 5A and 5B show examples of sequencing systems.

FIG. 6 shows an example of a computing device that can be used in, or inconjunction with, a sequencing machine and/or a computer system.

FIG. 7 shows the general architecture of the MT algorithm.

FIG. 8 shows pairwise analysis of nearby heterozygous SNPs.

FIG. 9 shows an example of the selection of an hypothesis and theassignment of a score to the hypothesis.

FIG. 10 shows graph construction.

FIG. 11 shows graph optimization.

FIG. 12 shows contig alignment.

FIG. 13 shows parent-assisted universal phasing.

FIG. 14 shows natural contig separations.

FIG. 15 shows universal phasing. Sequence legend: Mom: C-CGCAG TAGCTTACGAATCG (SEQ ID NO:1); Dad: G-ATTTA ACTGAGC ACTTGGC (SEQ ID NO:2).

FIG. 16 shows error detection using MT.

FIG. 17 shows an example of a method of decreasing the number of falsenegatives in which a confident heterozygous SNP call could be madedespite a small number of reads.

DETAILED DESCRIPTION

As used herein and in the appended claims, the singular forms “a,” “an,”and “the” include plural referents unless the context clearly dictatesotherwise. Thus, for example, reference to “a polymerase” refers to oneagent or mixtures of such agents, and reference to “the method” includesreference to equivalent steps and/or methods known to those skilled inthe art, and so forth.

Unless defined otherwise, all technical and scientific terms used hereinhave the same meaning as commonly understood by one of ordinary skill inthe art to which this invention belongs. All publications mentionedherein are incorporated herein by reference for the purpose ofdescribing and disclosing devices, compositions, formulations andmethodologies which are described in the publication and which might beused in connection with the presently described invention.

Where a range of values is provided, it is understood that eachintervening value, to the tenth of the unit of the lower limit unlessthe context clearly dictates otherwise, between the upper and lowerlimit of that range and any other stated or intervening value in thatstated range is encompassed within the invention. The upper and lowerlimits of these smaller ranges may independently be included in thesmaller ranges is also encompassed within the invention, subject to anyspecifically excluded limit in the stated range. Where the stated rangeincludes one or both of the limits, ranges excluding either both ofthose included limits are also included in the invention.

In the following description, numerous specific details are set forth toprovide a more thorough understanding of the present invention. However,it will be apparent to one of skill in the art that the presentinvention may be practiced without one or more of these specificdetails. In other instances, well-known features and procedures wellknown to those skilled in the art have not been described in order toavoid obscuring the invention.

Although the present invention is described primarily with reference tospecific embodiments, it is also envisioned that other embodiments willbecome apparent to those skilled in the art upon reading the presentdisclosure, and it is intended that such embodiments be contained withinthe present inventive methods.

The practice of the present invention may employ, unless otherwiseindicated, conventional techniques and descriptions of organicchemistry, polymer technology, molecular biology (including recombinanttechniques), cell biology, biochemistry, and immunology, which arewithin the skill of the art. Such conventional techniques includepolymer array synthesis, hybridization, ligation, and detection ofhybridization using a label. Specific illustrations of suitabletechniques can be had by reference to the example herein below. However,other equivalent conventional procedures can, of course, also be used.Such conventional techniques and descriptions can be found in standardlaboratory manuals such as Genome Analysis: A Laboratory Manual Series(Vols. I-IV), Using Antibodies: A Laboratory Manual, Cells: A LaboratoryManual, PCR Primer: A Laboratory Manual, and Molecular Cloning: ALaboratory Manual (all from Cold Spring Harbor Laboratory Press),Stryer, L. (1995) Biochemistry (4th Ed.) Freeman, New York, Gait,“Oligonucleotide Synthesis: A Practical Approach” 1984, IRL Press,London, Nelson and Cox (2000), Lehninger, Principles of Biochemistry 3rdEd., W. H. Freeman Pub., New York, N.Y. and Berg et al. (2002)Biochemistry, 5th Ed., W. H. Freeman Pub., New York, N.Y., all of whichare herein incorporated in their entirety by reference for all purposes.

Methods for Sequencing Complex Nucleic Acids—Overview

According to one aspect of the invention, methods are provided formultiple tagging of individual long fragments of target nucleic acids,or polynucleotides, including without limitation complex nucleic acids.Long fragments of a target nucleic acid or polynucleotide (e.g., 100 kbor longer) are tagged by a method that introduces a tag or barcode intomultiple sites in each fragment. Ideally, each fragment has introducedinto it multiple copies of one unique tag—a fragment-specific tag—or aunique pattern of insertion of multiple tags—a fragment-specific tagpattern. However, this is not necessarily the case: in some embodiments,a long fragment may have inserted into it more than one distinct tag,and more than one fragment may have inserted into them the same tag.After tagging, subfragments of the long fragments are produced, eachsubfragment including a tag. Commonly, the subfragments are amplified(e.g., by PCR). The subfragments, including the tags that are part ofeach subfragment, are then sequenced. The tag sequence permits thesequence data obtained from each subfragment to be assigned to the longfragment from which each subfragment is derived, which facilitatessequence mapping and assembly and the ordering of alleles (or hets) intoa haplotype of the target nucleic acids.

For accurate clinical sequencing and haplotyping of individual humangenomes from a small number of cells, long genomic fragments (˜100 kb orlonger) are preferable, although shorter fragments may be used. Assuming100 kb fragments, a human genome would have about 6×10⁴ fragments percell, and ˜18 cells would generate about 1 million fragments. DNA tagsthat are 12 bases long (12-mers) or longer have enough sequencediversity (16 million to over one billion) to tag each fragment with aunique tag.

We provide several methods for associating copies of the same long tagto hundreds of ˜1 kb sub-regions of ˜100 kb gnomic fragments in ahomogenous reaction without any physical compartmentalization (e.g.,droplets in an emulsion).

A description of various methods for practicing MT is provided below.

In some embodiments of the invention, such methods lead to a most (e.g.,50%, 60%, 70%, 80%, 90% or more) of the long fragments of a targetnucleic acid being tagged with multiple tag-containing sequences thatinclude the same tag sequence. Such methods minimize tagging withdifferent tag sequences, for example: selecting the proper ratio oftag-containing sequences to long fragments; selecting the properdilution or DNA concentration; minimizing molecule movement afterinitiation of the tagging process, for example, by mixing DNA fragments,tag-containing sequences, and enzymes and buffers at low temperature,waiting for liquid movements to stop, and then increasing thetemperature of the mixture in order to activate enzymatic processes);tethering a single tag-containing sequence to a single long DNA fragmentby covalent or non-covalent binding; and other techniques. There areseveral ways to attach or tether a single bead or nanoball with multiplecopies of a particular tag-containing sequence to a single long fragmentof a target nucleic acid. For example, a homopolymer sequence (e.g., andA-tail) may be added to the long fragment using terminal transferase oran adaptor with a selected sequence may be ligated to an end or ends ofthe long fragment. A complementary sequence may be added to the end ofor included within the tag-containing sequence or nanoball such that,under selected appropriate conditions, the tag-containing sequence ortag assembly anneals with the corresponding complementary sequence onthe long fragment. Preferably, a long fragment can anneal to only onetag-containing sequence or tag assembly.

In other embodiments of the invention, rather than maximizing the numberof long fragments with a single tag sequence inserted at multiplelocations along the long fragment, provide conditions under whichmultiple tags with different sequences are inserted at multiplelocations, creating a unique pattern or “fingerprint” for each longfragment that is provided by a unique pattern of insertion of thedifferent-sequence tags.

By sequencing 50% of each ˜1 kb fragment, ˜1× sequence coverage would begenerated for each genomic fragment because tagged fragments aregenerated from both strands of dsDNA. If we sequence 25% (½ readcoverage per fragment), we would observe the linkage of two regions in25% of fragments. For the same read budget we can increase the number offragments two-fold and have only two-fold reduction in the observedlinkages. For 25% read (125 bases form each end of 1 kb fragment) and 36starting cells we will observe nine linkages instead of ˜18 linkages for18 cells if we read 50% of DNA (250 bases from each end of a ˜1 kbfragment. If only ˜60 bases can be read from each fragment, it is betterto use 300-500 bp fragments that still make very useful mate-pairs.

If sequencing a fraction of DNA from each ˜1 kb fragment, we needexponentially more initial fragments. For example, if we sequenceone-half, we need 4× more fragments.

MT avoids subcloning of fragments of a complex nucleic acid into avector and subsequent replication in a host cell, or the need to isolateindividual chromosomes (e.g., metaphase chromosomes). It also does notrequire aliquoting fragments of a target nucleic acid. MT can be fullyautomated, making it suitable for high-throughput, cost-effectiveapplications. Tagging ˜1 kb sub-regions of long (˜100 kb or longer)genomic fragments with the same unique tag has many applications,including haplotyping diploid or polyploid genomes, efficient de novogenome sequence assembly, and error correction.

The advantages of MT include:

-   -   A practically unlimited number of individual DNA fragments can        be uniquely tagged, providing maximal information for de novo        assembly, for example.    -   MT may be performed in a single reaction vessel (e.g., tube,        well in a multi-well plate, etc.) in a small number of steps and        is easy to scale and automate; there is no need for a large        number of aliquots or nanodrops.    -   DNA amplification of shorter DNA fragments (e.g., tagged        subfragments) by PCR results in less coverage bias and no loss        of sequences from fragment ends compared to multiple        displacement amplification (MDA) of long fragments, for example.    -   One MT method, which employs nicking and a primer-ligation        process, uses both strand of the dsDNA, thereby doubling the        sequence coverage per fragment (longer mate-pairs for the same        read-length).    -   Reduce dramatically the computational demands and associated        costs of sequence mapping and assembly.    -   Substantial reduction in errors or questionable base calls that        can result from current sequencing technologies, including, for        example, systematic errors that are characteristic of a given        sequencing platform or mutations introduced by DNA        amplification. MT thereby provides a highly accurate sequence of        a human genome or other complex nucleic acid, minimizing the        need for follow up confirmation of detected variants and        facilitates adoption of human genome sequencing for diagnostic        applications.

MT can be used as a preprocessing method with any known sequencingtechnology, including both short-read and longer-read methods. MT alsocan be used in conjunction with various types of analysis, including,for example, analysis of the transcriptome, methylome, etc. Because itrequires very little input DNA, MT can be used for sequencing andhaplotyping one or a small number of cells, which can be particularlyimportant for cancer, prenatal diagnostics, and personalized medicine.This can facilitate the identification of familial genetic disease, etc.By making it possible to distinguish calls from the two sets ofchromosomes in a diploid sample, MT also allows higher confidencecalling of variant and non-variant positions at low coverage. Additionalapplications of MT include resolution of extensive rearrangements incancer genomes and full-length sequencing of alternatively splicedtranscripts.

MT can be used to process and analyze complex nucleic acids, includingbut not limited to genomic DNA, that is purified or unpurified,including cells and tissues that are gently disrupted to release suchcomplex nucleic acids without shearing and overly fragmenting suchcomplex nucleic acids.

In one aspect, MT produces virtual read lengths of approximately100-1000 kb or longer in length, for example.

In addition to being applicable to all sequencing platforms, MT-basedsequencing is suitable for a wide variety of applications, includingwithout limitation the study of structural rearrangements in cancergenomes, full methylome analysis including the haplotypes of methylatedsites, and de novo assembly applications for metagenomics or novelgenome sequencing, even of complex polyploid genomes like those found inplants.

MT provides the ability to obtain actual sequences of individualchromosomes as opposed to just the consensus sequences of parental orrelated chromosomes (in spite of their high similarities and presence oflong repeats and segmental duplications). To generate this type of data,the continuity of sequence is in general established over long DNAranges.

A further aspect of the invention includes software and algorithms forefficiently utilizing MT data for whole chromosome haplotype andstructural variation mapping and false positive/negative errorcorrection.

In a further aspect, MT techniques of the invention reduce thecomplexity of DNA to be sequenced. Complexity reduction and haplotypeseparation in >100 kb long DNA can be helpful in more efficiently andcost effective sequence assembly and detection of sequence variations inhuman and other diploid and polyploid genomes.

Preparing Long Nucleic Acid Fragments

Target nucleic acids, including but not limited to complex nucleicacids, are isolated using conventional techniques, for example asdisclosed in Sambrook and Russell, Molecular Cloning: A LaboratoryManual, cited supra. In some cases, particularly if small amounts of thenucleic acids are employed in a particular step, it is advantageous toprovide carrier DNA, e.g., unrelated circular synthetic double-strandedDNA, to be mixed and used with the sample nucleic acids whenever onlysmall amounts of sample nucleic acids are available and there is dangerof losses through nonspecific binding, e.g., to container walls and thelike.

According to some embodiments of the invention, genomic DNA or othercomplex nucleic acids are obtained from an individual cell or smallnumber of cells with or without purification, by any known method.

Long fragments are desirable for the methods of the present invention.Long fragments of genomic DNA can be isolated from a cell by any knownmethod. A protocol for isolation of long genomic DNA fragments fromhuman cells is described, for example, in Peters et al., Nature487:190-195 (2012). In one embodiment, cells are lysed and the intactnuclei are pelleted with a gentle centrifugation step. The genomic DNAis then released through proteinase K and RNase digestion for severalhours. The material can be treated to lower the concentration ofremaining cellular waste, e.g., by dialysis for a period of time (i.e.,from 2-16 hours) and/or dilution. Since such methods need not employmany disruptive processes (such as ethanol precipitation,centrifugation, and vortexing), the genomic nucleic acid remains largelyintact, yielding a majority of fragments that have lengths in excess of150 kilobases. In some embodiments, the fragments are from about 5 toabout 750 kilobases in lengths. In further embodiments, the fragmentsare from about 150 to about 600, about 200 to about 500, about 250 toabout 400, and about 300 to about 350 kilobases in length. The smallestfragment that can be used for haplotyping is one containing at least twohets (approximately 2-5 kb); there is no maximum theoretical size,although fragment length can be limited by shearing resulting frommanipulation of the starting nucleic acid preparation.

Once the DNA is isolated, it is advantageous to avoid loss of sequencesfrom the ends of each fragment, since loss of such material can resultin gaps in the final genome assembly. In one embodiment, sequence lossis avoided through use of an infrequent nicking enzyme, which createsstarting sites for a polymerase, such as phi29 polymerase, at distancesof approximately 100 kb from each other. As the polymerase creates a newDNA strand, it displaces the old strand, creating overlapping sequencesnear the sites of polymerase initiation. As a result, there are very fewdeletions of sequence.

A controlled use of a 5′ exonuclease (either before or duringamplification) can promote multiple replications of the original DNAfrom a single cell and thus minimize propagation of early errors throughcopying of copies.

In other embodiments, long DNA fragments are isolated and manipulated ina manner that minimizes shearing or absorption of the DNA to a vessel,including, for example, isolating cells in agarose in agarose gel plugs,or oil, or using specially coated tubes and plates.

Fragmented DNA from a single cell can be duplicated by ligating anadaptor with single stranded priming overhang and using anadaptor-specific primer and phi29 polymerase to make two copies fromeach long fragment. This can generate four cells-worth of DNA from asingle cell.

Preparation of Tags

A “barcode” or “tag” is generally a unique sequence of nucleotides.Tagging long fragments as described involves introducing into longfragments multiple copies of sequences (adaptors, transposons, etc.)that include tags. Such introduced sequences are spaced apart on thefragment; the average spacing between adjacent introduced sequences isselected to permit the creation of subfragments of the long fragments byany preferred method, e.g.: by PCR amplification using primers that haveprimer binding sites in adjacent introduced sequences; by restrictiondigestion; or by other methods known in the art. Subsequently, sequencereads are generated by sequencing subfragments of the tagged longfragments. Such sequence reads can be assigned to the individual longfragment from which they are ultimately derived.

According to one embodiment, the average spacing between adjacentintroduced sequences is 100 bp, 200 bp, 300 bp, 400 bp, 500 bp, 600 bp,700 bp, 800 bp, 900 bp, 1000 bp, 1500 bp, 2000 bp, 2500 bp, 3000 bp,3500 bp, 4000 bp, or 5000 bp. According to another embodiment, theaverage spacing is between about 100 bp and about 5000 bp, or betweenabout 200 bp and about 4000 bp, or between about 300 bp and about 3000bp, or between about 300 bp and about 2000 bp, or between about 300 bpand about 1000 bp.

According to one embodiment, a majority of the subfragments, or 60%,70%, 80%, 90% or more, or substantially all of the subfragments, includea tag sequence.

According to one embodiment, a barcode- or tag-containing sequence isused that has two, three or more segments, one or more regions of knownsequence and one or more regions of degenerate sequence that serves asthe barcode(s) or tag(s). The known sequence (B) may include, forexample, PCR primer binding sites, transposon ends, restrictionendonuclease recognition sequences (e.g., sites for rare cutters, e.g.,Not I, Sac II, Mlu I, BssH II, etc.), or other sequences. The degeneratesequence (N) that serves as the tag is long enough to provide apopulation of different-sequence tags that is equal to or, preferably,greater than, the number of fragments of a target nucleic acid to beanalyzed.

According to one embodiment, the tag-containing sequence comprises oneregion of known sequence of any selected length (including, withoutlimitation. According to another embodiment the tag-containing sequencecomprises two regions of known sequence of a selected length that flanka region of degenerate sequence of a selected length, i.e.,B_(n)N_(n)B_(n), where N may have any length sufficient for tagging longfragments of a target nucleic acid, including, without limitation, N=10,11, 12, 13, 14, 15, 16, 17, 18, 19 or 20, and B may have any length thataccommodates desired sequences such as transposon ends, primer bindingsites, etc. For example, such an embodiment may be B₂₀N₁₅B₂₀.

Tag Synthesis on Beads.

A population of tag-containing DNA sequences can be PCR amplified onbeads in an water-in-oil (w/o) emulsion by known methods. See, e.g.,Tawfik and Griffiths Nature Biotechnology 16: 652-656 (1998); Dressmanet al., Proc. Natl. Acad. Sci. USA 100:8817-8820, 2003; and Shendure etal., Science 309:1728-1732 (2005). This results in many copies of eachsingle tag-containing sequence on each bead.

Concatemers (DNA Nanoballs) of Tag Sequences.

A circular or circularized (e.g., as with padlock probes) DNA templatecan be amplified by rolling circle replication (RCR). RCR uses the phi29 DNA polymerase, which is highly processive. The newly synthesizedstrand is released from the circular template, resulting in a longsingle-stranded DNA concatemer comprising many head-to-tail copies ofthe circular DNA template. The concatemer folds into a substantiallyglobular ball of DNA that is called a DNA nanoball (DNB). The length ofthe DNB and the number of copies of the DNA template can be controlledby the length of the RCR reaction. The nanoballs remain separated fromeach other in solution.

In one embodiment, a two- or three-segment design is utilized for thebarcodes used to tag long fragments. This design allows for a widerrange of possible barcodes by allowing combinatorial barcode segments tobe generated by ligating different barcode segments together to form thefull barcode segment or by using a segment as a reagent inoligonucleotide synthesis. This combinatorial design provides a largerrepertoire of possible barcodes while reducing the number of full-sizebarcodes that need to be generated. In further embodiments, uniqueidentification of each long fragment is achieved with 8-12 base pair (orlonger) barcodes.

In one embodiment, two different barcode segments are used. A and Bsegments are easily be modified to each contain a different half-barcodesequence to yield thousands of combinations. In a further embodiment,the barcode sequences are incorporated on the same adapter. This can beachieved by breaking the B adaptor into two parts, each with a halfbarcode sequence separated by a common overlapping sequence used forligation. The two tag components have 4-6 bases each. An 8-base (2×4bases) tag set is capable of uniquely tagging 65,000 sequences. Both 2×5base and 2×6 base tags may include use of degenerate bases (i.e.,“wild-cards”) to achieve optimal decoding efficiency.

In further embodiments, unique identification of each sequence isachieved with 8-12 base pair error correcting barcodes.

According to one embodiment of the invention, one starts with more longfragments than are needed for sequencing to achieve adequate sequencecoverage and tags only a only a portion of the long fragments with alimited number of tag-containing sequences, or tag assemblies—whichinclude many, perhaps hundreds, of copies of one tag sequence—toincrease the probability of unique tagging of the long fragments.Non-tagged subfragments lack introduced sequences that provideprimer-binding or capture-oligo binding and may be eliminated indownstream processing. Such tag assemblies include, for example,end-to-end concatemers of tag-containing sequences created by rollingcircle replication (DNA nanoballs), beads to which are attached manycopies of the tag-containing sequences, or other embodiments.

According to another embodiment, in order to obtain uniform genomecoverage in the case of samples with a small number of cells (e.g., 1,2, 3, 4, 5, 10, 10, 15, 20, 30, 40, 50 or 100 cells from a microbiopsyor circulating tumor or fetal cells, for example), all long fragmentsobtained from the cells are tagged.

Tagging Single Long Fragments

The methods of the present invention employ various approaches tointroduce multiple copies of a tag at multiple spaced-apart sites alonga long fragment (e.g., 100 kb or longer) of the target nucleic acidwithout the need to divide the long fragments into aliquots (as in thelong fragment read technology): the entire process can be performed in asingle tube or well in a microtiter plate.

According to one embodiment of the invention, tags are introduced atintervals of between about 300 bp and 1000 bp along the fragment. Thisspacing can be shorter or longer, depending on the desired fragment sizefor subsequent processing, e.g., library construction and sequencing.After tagging, each subfragment of the long fragment and any sequenceinformation derived from it can be assigned to a single long fragment.

Various exemplary methods for tagging single fragments are shown inFIGS. 1-4.

(1) Tagging with Transposons.

As shown in FIG. 1, a first approach involves in vitro transposition. Apopulation of tagged transposons is used, i.e., DNA constructs thatinclude transposon ends and, near each of the ends, unique tag sequences(the same tag sequence near both ends) and a common PCR primer bindingsite. The population of transposons is combined with long fragments of atarget nucleic acid. Addition of transposase causes in vitrotransposition of several of the tagged transposons into the longfragments. Each long fragment has a unique pattern of transposoninsertion, and each inserted transposon has a unique bar code. Inaddition, the act of transposition replicates 9 bp of sequence at eachend of the transposon that further distinguishes each transposoninsertion event (and may be considered another form of “tagging”).

PCR is performed using primers that bind to the PCR primer binding sitesof each inserted transposon. The resulting PCR amplicons include aportion of the long fragment that lies between adjacent transposons; ateach end such amplicons include sequences from the end of an adjacenttransposon, including the unique tag sequence for that transposon. Aftersequencing the PCR amplicons, it is possible not only to map thesequence reads to a reference genome, assuming such is available, but touse the tags to build contigs to guide de novo assembly, as shown inFIG. 1. Each sequence read is associated with a tag sequence, and thetag sequence (or pattern of tags) corresponds to a single fragment.Thus, sequence reads from the same fragment should map within the sameregion of the target nucleic acid. In most cases two different ampliconshave the same unique tag from one transposon at their ends and thusadjoin one another in the large fragment from which they are derived. Ifmore than one genome equivalent of long fragments is analyzed (e.g., 2,3, 4, 5, 10, or 20 or more genome equivalents) building up contigs outof sequence reads derived from overlapping long fragments isstraightforward.

(2) Tagging with Hairpins.

As shown in FIGS. 2A and 2B, this approach begins with long fragments ofa target nucleic acid that are denatured to form two complementarysingle strands from each fragment. It also uses a population ofoligonucleotides that form hairpins, each including tag sequences in theloop that flank PCR primer binding sites and having a short stretch ofrandom bases (3-5 bases) at each end. The hairpin oligos are annealed tothe single stranded form of the starting long fragments spaced apart,for example, by about 300 to 1000 bp. Each long fragment has a uniquepattern of annealed hairpins. After annealing the single stranded regionbetween adjacent hairpins is filled in with a 5′-3′ polymerase thatlacks strand displacement, followed by ligase treatment to seal theremaining nick. PCR amplification using primers that bind to the PCRbinding sites between the bar code sequences of each hairpin createsamplicons that have a portion of the long fragment that lies betweenadjacent hairpin oligonucleotides; at each end such amplicons includesequences from the loop of an adjacent hairpin oligonucleotide,including the unique tag sequence for that transposon. As for method (1)above, the bar code sequences at the ends of the PCR amplicons can beused to build contigs to guide de novo mapping and assembly.

(3) Tagging with Transposons on a DNA Nanoball or Bead.

As shown in FIG. 3, this approach employs beads covered with transposonsequences or a concatemer of transposon sequences created by rollingcircle replication of a circular DNA that includes the transposonsequence—a transposon nanoball. As in method (1) above, the “transposonsequences” are DNA constructs that include (i) transposon ends and, at aselected location in the transposon sequences between the transposonends, e.g., near each of the transposon ends, (ii) unique tag sequences(the same tag sequence near both ends), and (iii) a common PCR primerbinding site. The transposon-containing bead or nanoball is combinedwith the long fragments of a target nucleic acid. Conditions areselected to promote the interaction of only a tag assembly, i.e., beador nanoball bearing a single transposon sequence, with each longfragment. For example, at the correct dilution, only one bead ornanoball interacts with each long fragment in most cases, sincediffusion is slow and most transposons don't travel far from a longfragment. Alternatively, the transposon sequence or another sequence onthe transposon assembly (e.g., an adaptor ligated to an end of thetransposon sequence or concatemer; a homopolymer sequence added by aterminal transferase) Upon addition of transposase, transpositionoccurs. In most cases, each fragment has multiple copies of the sametransposon. A minority of the long fragments receive copies of more thanone transposon. Also, in some cases, a transposon with a particular tagtransposes into more than one long fragment. As in method (1), PCRamplification is performed using primers that bind to the PCR primerbinding sites of each inserted transposon. The resulting PCR amplicons(between about 300 bp and 1000 bp in length) include a portion of thelong fragment that lies between adjacent transposons; at each end suchamplicons include sequences from the end of an adjacent transposon,including the unique tag sequence for that transposon. After sequencing,sequence reads are mapped and assembled.

In this method, because most long fragments are tagged with multiplecopies of a single transposon, the resulting amplicons have the same tagat each end. The tags permit each sequence read to be associated withthe same long fragment, although it is not possible to build up contigsbased on the ordering of the tag sequences alone as in methods (1) and(2). If more than one transposon inserts into a single long fragment, itis most likely that all of the transposons that insert into one longfragment insert only into that one long fragment and not into otherfragments. As a result, sequence reads associated with each of theinserted tags maps closely together in the genome (or other targetnucleic acid). Even if this is not the case, and the same transposonjumps into more than one fragment, the likelihood is high that thefragments into which such transposon is inserted are non-overlapping, inwhich case the resulting sequence reads map to widely separated regionsof the genome. Mapping and assembly software can account for theseevents and correctly map and assemble the sequence reads into a genomesequence and order sequence polymorphisms (hets) into a haplotype.

Tagging with an Tagged Adaptor Carried on a DNA Nanoball or Bead.

In this method, shown in FIG. 4A, long double-stranded fragments of agenome (or other target nu acid) are nicked at random locations on bothstrands using an agent such as DNase I that nicks DNA double strands(i.e., a “nickase”) and DNA polymerase I large (Klenow) fragment, whichretains polymerization and 3′→5′ exonuclease activity, but has lost5′→3′ exonuclease activity. No dNTPs are included in the reaction. A 3′common adaptor is ligated to the 3′ end of each strand at a nick. A beador nanoball with many copies of a sequence that includes (i) the tagsequence and (ii) a sequence that is complementary to the 3′ commonadapter is added under conditions that permit the 3′ common adaptor tohybridize to the complementary sequence. At the proper ratio of longfragments to beads or nanoballs and at the proper dilution, most of thefragments are spatially associated with one (or less frequently 2 ormore) beads or nanoballs, and copies of the 3′ common adaptor hybridizeto the complementary sequence on a single bead or nanoball, since asingle hybridization event leads to a physical interaction between thelong fragment and the bead or nanoball, bringing other complementarysequences into close proximity. Alternatively, one can use complementaryDNA sequences such as an A-tail on the long fragment and a T-tail orpoly-T region on the tag-containing sequences, or other interactivemoieties, to force the interaction of the long fragment andtag-containing sequences in order to increase the likelihood that eachlong fragment has introduced into it multiple copies of a singletag-containing sequence. Next, the tag-containing nucleic acids on thebead or nanoball are fragmented, e.g., with a restriction endonuclease,which results in common adaptors ligated to the long fragmenthybridizing to complementary sequences that are included in the nucleicacids released from the bead or nanoball. Primer extension using DNApolymerase I large fragment (Klenow) or a similar DNA polymerase resultsin the creation of a 3′ tagged molecule spaced apart on the longfragment every 300-1000 bp. The long DNA molecule can now be denaturedand an oligonucleotide can be hybridized to the 3′ common adaptor;extension with Klenow fragment or a similar polymerase results in ablunt-ended, double-stranded DNA molecule that can be ligated to a 5′common adaptor and PCR amplified. The resulting PCR amplicons(effectively tagged subfragments of the long DNA fragments) are thensequenced, mapped and assembled in a fashion similar to that describedin method (3).

Thus, according to this method of the invention, the MT processcomprises:

(1) “Clonal” copying of barcode templates and required adapters, forexample by (a) rolling circle replication (RCR) to make a concatemerwith hundreds of copies of the same tag or by (b) emulsion PCR on beadsto create thousands of copies. Optionally, the copied unit may representa transposon.

(2) Mixing long genomic fragments and tag-adapter concatemers or beadsin the proper ratio and in proper concentrations to have majority, mostor almost all genomic fragments spatially associated with one concatemerand infrequently with two or more.

(3) Adding a universal primer to genomic DNA by: (a) nicking of genomicDNA at predefined frequency (e.g., 1 kb) using partial nicking withfrequent nicker or other methods; controlled nick translation can beused to further randomize fragment start sites; optionally a small gapmay be created at the nicking site, e.g., by exo activity of Pol I orKlenow without dNTPs; (b) ligating a primer by 5′ end to 3′ end ofnicked DNA by providing the primer hybridized with a short complementarydideoxy oligo at the 5′ end; this primer is complementary to an adapternext to the barcode. Optionally this step can be done before step two ormixing genomic DNA with clonal tags;

(4) Copying the tag from tag donor (DNA nanoball or bead) and anotheradapter by primer extension using tag templates. After DNA denaturingthis results in ˜1 kb ssDNA fragments with an adapter-barcode-adapterextension at 3′ end. These fragments can be used as sequencing templatesby a primer complementary to the 3′ end adapter or converted in dsDNA bythe same primer and further process (e.g., ligate an adapter on theother end, amplify, circularize) before sequencing.

Optionally steps 3 and 4 can be replaced by transposon insertion andfragmenting or amplification if concatemers or beads represent clones oftagged transposons.

FIG. 4B shows an alternative approach to inserting tag-containingsequences that does not rely on nicking. Long fragments are denatured(e.g., by heating) to produce complementary single strands. Randomprimers (N-mers) are annealed to the single strands and extended withpolymerase. An alkaline phosphatase (e.g., shrimp alkaline phosphatase,SAP) is added, and polymerase having a 3′→5′ exonuclease function (e.g.,Klenow) is used to create gaps, and the resulting partiallydouble-stranded product is handled as described above and in FIG. 4A,beginning with 3′ ligation of a common adaptor.

FIG. 4C-D show a second alternative approach to inserting tag-containingsequences that does not rely on nicking. In this approach, twooligonucleotides are annealed to a tag-containing sequence (carried on abead or, as shown, as a monomer unit of a DNA concatemer or nanoball):(i) a common primer, which is annealed upstream of the tag or barcodesequence, and (ii) a common adapter that is annealed downstream of thetag. The primer is extended and ligase is added to ligate the primerextension product to the common adaptor. This ligation product thusincludes the tag sequence and, at its 3′ end, the common adaptor. Anpopulation of oligonucleotides that includes (i) a degenerate sequence(random N-mer) at its 5′ end, (ii) a sequence complementary to thecommon adaptor, and (iii) noncomplementary sequence (not shown in FIG.4C) is annealed to the ligation product from the previous step and aprimer extension is performed, adding to the 3′ end of the ligationproduct a degenerate sequence complementary to that on eacholigonucleotide (which is subsequently removed, e.g., chewed away). Theresulting product (a population of “tagged adaptors” each with adegenerate sequence at their 3′ ends) is then released from the bead ornanoball, e.g., by heat denaturation. The tagged adaptors are annealedto a single strand of the long fragment (produced by denaturing thedouble stranded long fragment); as shown in FIG. 4D, the differentdegenerate sequences at the ends of various tagged adaptors anneal tocomplementary sequences spaced apart along the long fragment. Asdescribed above, a polymerase is added to extend the tagged adaptor, andthe extension product includes a sequence complementary to a region ofthe long fragment. The resulting molecules, which include a taggedadaptor joined to a sequence from the long fragment can then be used tocreate tagged subfragments of the long fragment as described above (FIG.4A).

The method described in FIGS. 4A and 4B results in PCR amplicons thatare, effectively, tagged subfragments of the long DNA fragments. This isadvantageous if short-read sequencing methods are used. There are avariety of ways to create such a series of fragments.

For example, it is possible to create a series of such subfragments withshorter and shorter regions of the long DNA fragments as shown in FIG.4E. Starting with the blunt-ended primer-extended tagged subfragmentthat results from PCR amplification, a 3′ adaptor is ligated to thetagged subfragments. One end of the adaptor includes an overhang; theother end is a blunt end that includes a blocked nucleotide (e.g., addNTP). After ligation of the 3′ adaptor, the subfragment is denaturedand another round of primer extension is performed using controlled nicktranslation. The primer extension is stopped before completion such thatthe primer does not extend all the way to the end of the complementarystrand. A 3′ adaptor is ligated to the end of the extended strand. Thisprocess can be repeated as many times as desired with the extent ofprimer extension varied in order to create a series of fragments havinga common 5′ end that are shortened on their 3′ ends. Details of theblocked adaptor strategy and of controlled nick translation areprovided, for example, in U.S. patent application Ser. No. 12/329,365(published as U.S. 2012-0100534 A1) and Ser. No. 12/573,697 (publishedas US-2010-0105052-A1).

Another approach to creating a series of such subfragments with shorterand shorter regions of the long DNA fragments as shown in FIG. 4F. Thisapproach also uses controlled nick translation. Subfragments arecircularized then split into two or more separate wells. Controlled nicktranslation is performed to a different extent in the various wells inorder to create subfragments with a common 5′ end that are shortened ontheir 3′ ends to various degrees. The subfragments can then be pooledand the process continued.

Producing Subfragments of Tagged Long Fragments

After tagging, the long fragments of the target nucleic acid aresubfragmented to a desired size by amplification (e.g., by PCR),restriction enzyme digestion (e.g., using a rare cutter that has arecognition site within a tag-containing sequence introduced into longfragments), or by other conventional techniques, including enzymaticdigestion, shearing, sonication, etc.

Subfragment sizes can vary depending on the source target nucleic acidand the library construction methods used, but for standard whole-genomesequencing such fragments typically range from 50 to 2000 nucleotides inlength. In another embodiment, the fragments are 300 to 600 or 200 to2000 nucleotides in length. In yet another embodiment, the fragments are10-100, 50-100, 50-300, 100-200, 200-300, 50-400, 100-400, 200-400,300-400, 400-500, 400-600, 500-600, 50-1000, 100-1000, 200-1000,300-1000, 400-1000, 500-1000, 600-1000, 700-1000, 700-900, 700-800,800-1000, 900-1000, 1500-2000, 1750-2000, and 50-2000 nucleotides inlength.

In a further embodiment, fragments of a particular size or in aparticular range of sizes are isolated. Such methods are well known inthe art. For example, gel fractionation can be used to produce apopulation of fragments of a particular size within a range ofbasepairs, for example for 500 base pairs+50 base pairs.

Starting with about 5 to about 1,000,000 genome-equivalents of longfragment DNA ensure that the population of long fragments covers theentire genome. Libraries containing nucleic acid templates generatedfrom such a population of overlapping fragments will provide most or allof the sequence of an entire genome.

Amplification

Before or after any step outlined herein, an amplification step can beused to ensure that enough of the nucleic acid is available forsubsequent steps. According to one embodiment of the invention, methodsare provided for sequencing small quantities of complex nucleic acids,including those of higher organisms, in which such complex nucleic acidsare amplified in order to produce sufficient nucleic acids forsequencing by the methods described herein. A single human cell includesapproximately 6.6 picograms (pg) of genomic DNA. Sequencing of complexnucleic acids of a higher organism can be accomplished using 1 pg, 5 pg,10 pg, 30 pg, 50 pg, 100 pg, or 1 ng or more of a complex nucleic acidas the starting material, which is amplified by any nucleic acidamplification method known in the art, to produce, for example, 200 ng,400 ng, 600 ng, 800 ng, 1 μg, 2 μg, 3 μg, 4 μg, 5 μg, 10 μg or greaterquantities of the complex nucleic acid. We also disclose nucleic acidamplification protocols that minimize GC bias. However, the need foramplification and subsequent GC bias can be reduced further simply byisolating one cell or a small number of cells, culturing them for asufficient time under suitable culture conditions known in the art, andusing progeny of the starting cell or cells for sequencing.

Such amplification methods include without limitation: multipledisplacement amplification (MDA), polymerase chain reaction (PCR),ligation chain reaction (sometimes referred to as oligonucleotide ligaseamplification OLA), cycling probe technology (CPT), strand displacementassay (SDA), transcription mediated amplification (TMA), nucleic acidsequence based amplification (NASBA), rolling circle amplification (RCA)(for circularized fragments), and invasive cleavage technology.

Amplification can be performed after fragmenting or before or after anystep outlined herein.

MDA amplification protocol with Reduced GC Bias.

In one aspect, the present invention provides methods of nucleic acidamplification in which the nucleic acid is faithfully amplified, e.g.,approximately 30,000-fold depending on the amount of starting DNA.

According to one embodiment of MT methods of the present invention, MTbegins with treatment of genomic nucleic acids, usually genomic DNA,with a 5′ exonuclease to create 3′ single-stranded overhangs. Suchsingle stranded overhangs serve as MDA initiation sites. Use of theexonuclease also eliminates the need for a heat or alkaline denaturationstep prior to amplification without introducing bias into the populationof fragments. In another embodiment, alkaline denaturation is combinedwith the 5′ exonuclease treatment, which results in a reduction in biasthat is greater than what is seen with either treatment alone. Thefragments are then amplified.

In one embodiment, a phi29-based multiple displacement amplification(MDA) is used. Numerous studies have examined the range of unwantedamplification biases, background product formation, and chimericartifacts introduced via phi29 based MDA, but many of these shortcomings have occurred under extreme conditions of amplification (greaterthan 1 million fold). Commonly, MT employs a substantially lower levelof amplification and starts with long DNA fragments (e.g., ˜100 kb),resulting in efficient MDA and a more acceptable level of amplificationbiases and other amplification-related problems.

We have developed an improved MDA protocol to overcome problemsassociated with MDA that uses various additives (e.g., DNA modifyingenzymes, sugars, and/or chemicals like DMSO), and/or differentcomponents of the reaction conditions for MDA are reduced, increased orsubstituted to further improve the protocol. To minimize chimeras,reagents can also be included to reduce the availability of thedisplaced single stranded DNA from acting as an incorrect template forthe extending DNA strand, which is a common mechanism for chimeraformation. A major source of coverage bias introduced by MDA is causedby differences in amplification between GC-rich verses AT-rich regions.This can be corrected by using different reagents in the MDA reactionand/or by adjusting the primer concentration to create an environmentfor even priming across all % GC regions of the genome. In someembodiments, random hexamers are used in priming MDA. In otherembodiments, other primer designs are utilized to reduce bias. Infurther embodiments, use of 5′ exonuclease before or during MDA can helpinitiate low-bias successful priming, particularly with longer (i.e.,200 kb to 1 Mb) fragments that are useful for sequencing regionscharacterized by long segmental duplication (i.e., in some cancer cells)and complex repeats.

In some embodiments, improved, more efficient fragmentation and ligationsteps are used that reduce the number of rounds of MDA amplificationrequired for preparing samples by as much as 10,000 fold, which furtherreduces bias and chimera formation resulting from MDA.

In some embodiments, the MDA reaction is designed to introduce uracilsinto the amplification products in preparation for CoRE fragmentation.In some embodiments, a standard MDA reaction utilizing random hexamersis used to amplify the fragments in each well; alternatively, random8-mer primers can be used to reduce amplification bias (e.g., GC-bias)in the population of fragments. In further embodiments, severaldifferent enzymes can also be added to the MDA reaction to reduce thebias of the amplification. For example, low concentrations ofnon-processive 5′ exonucleases and/or single-stranded binding proteinscan be used to create binding sites for the 8-mers. Chemical agents suchas betaine, DMSO, and trehalose can also be used to reduce bias.

After amplification of the nucleic acids in a sample, the amplificationproducts may optionally be fragmentated. In some embodiments the CoREmethod is used to further fragment the fragments followingamplification. In such embodiments, MDA amplification of fragments isdesigned to incorporate uracils into the MDA products. The MDA productis then treated with a mix of Uracil DNA glycosylase (UDG), DNAglycosylase-lyase Endonuclease VIII, and T4 polynucleotide kinase toexcise the uracil bases and create single base gaps with functional 5′phosphate and 3′ hydroxyl groups. Nick translation through use of apolymerase such as Taq polymerase results in double stranded blunt-endbreaks, resulting in ligatable fragments of a size range dependent onthe concentration of dUTP added in the MDA reaction. In someembodiments, the CoRE method used involves removing uracils bypolymerization and strand displacement by phi29. The fragmenting of theMDA products can also be achieved via sonication or enzymatic treatment.Enzymatic treatment that could be used in this embodiment includeswithout limitation DNase I, T7 endonuclease I, micrococcal nuclease, andthe like.

Following fragmentation of the MDA products, the ends of the resultantfragments may be repaired. Many fragmentation techniques can result intermini with overhanging ends and termini with functional groups thatare not useful in later ligation reactions, such as 3′ and 5′ hydroxylgroups and/or 3′ and 5′ phosphate groups. It may be useful to havefragments that are repaired to have blunt ends. It may also be desirableto modify the termini to add or remove phosphate and hydroxyl groups toprevent “polymerization” of the target sequences. For example, aphosphatase can be used to eliminate phosphate groups, such that allends contain hydroxyl groups. Each end can then be selectively alteredto allow ligation between the desired components. One end of thefragments can then be “activated” by treatment with alkalinephosphatase.

Nucleic Acid Sequencing

MT methods described herein can be used as a pre-processing step forsequencing diploid genomes using any sequencing method known in the art,including for example without limitation, polymerase-basedsequencing-by-synthesis (e.g., HiSeq 2500 system, Illumina, San Diego,Calif.), ligation-based sequencing (e.g., SOLiD 5500, Life TechnologiesCorporation, Carlsbad, Calif.), ion semiconductor sequencing (e.g., IonPGM or Ion Proton sequencers, Life Technologies Corporation, Carlsbad,Calif.), zero-mode waveguides (e.g., PacBio RS sequencer, PacificBiosciences, Menlo Park, Calif.), nanopore sequencing (e.g., OxfordNanopore Technologies Ltd., Oxford, United Kingdom), pyrosequencing(e.g., 454 Life Sciences, Branford, Conn.), or other sequencingtechnologies. Some of these sequencing technologies are short-readtechnologies, but others produce longer reads, e.g., the GS FLX+ (454Life Sciences; up to 1000 bp), PacBio RS (Pacific Biosciences;approximately 1000 bp) and nanopore sequencing (Oxford NanoporeTechnologies Ltd.; 100 kb). For haplotype phasing, longer reads areadvantageous, requiring much less computation, although they tend tohave a higher error rate and errors in such long reads may need to beidentified and corrected according to methods set forth herein beforehaplotype phasing.

According to one embodiment, sequencing is performed using combinatorialprobe-anchor ligation (cPAL) as described, for example, in U.S. PatentApplication Publications 2010/0105052 and US2007099208, and U.S. patentapplication Ser. Nos. 13/448,279, 13/447,087, 11/679,124 (published asUS 2009/0264299); Ser. No. 11/981,761 (US 2009/0155781); Ser. No.11/981,661 (US 2009/0005252); Ser. No. 11/981,605 (US 2009/0011943);Ser. No. 11/981,793 (US 2009-0118488); Ser. No. 11/451,691 (US2007/0099208); Ser. No. 11/981,607 (US 2008/0234136); Ser. No.11/981,767 (US 2009/0137404); Ser. No. 11/982,467 (US 2009/0137414);Ser. No. 11/451,692 (US 2007/0072208); Ser. No. 11/541,225 (US2010/0081128; Ser. No. 11/927,356 (US 2008/0318796); Ser. No. 10/927,388(US 2009/0143235); Ser. No. 11/938,096 (US 2008/0213771); Ser. No.11/938,106 (US 2008/0171331); Ser. No. 10/547,214 (US 2007/0037152);Ser. No. 11/981,730 (US 2009/0005259); Ser. No. 11/981,685 (US2009/0036316); Ser. No. 11/981,797 (US 2009/0011416); Ser. No.11/934,695 (US 2009/0075343); Ser. No. 11/934,697 (US 2009/0111705);Ser. No. 11/934,703 (US 2009/0111706); Ser. No. 12/265,593 (US2009/0203551); Ser. No. 11/938,213 (US 2009/0105961); Ser. No.11/938,221 (US 2008/0221832); Ser. No. 12/325,922 (US 2009/0318304);Ser. No. 12/252,280 (US 2009/0111115); Ser. No. 12/266,385 (US2009/0176652); Ser. No. 12/335,168 (US 2009/0311691); Ser. No.12/335,188 (US 2009/0176234); Ser. No. 12/361,507 (US 2009/0263802),Ser. No. 11/981,804 (US 2011/0004413); and Ser. No. 12/329,365;published international patent application numbers WO2007120208,WO2006073504, and WO2007133831, all of which are incorporated herein byreference in their entirety for all purposes.

Exemplary methods for calling variations in a polynucleotide sequencecompared to a reference polynucleotide sequence and for polynucleotidesequence assembly (or reassembly), for example, are provided in U.S.patent publication No. 2011-0004413, (application Ser. No. 12/770,089)which is incorporated herein by reference in its entirety for allpurposes. See also Drmanac et al., Science 327, 78-81, 2010. Alsoincorporated by references in their entirety and for all purposes arecopending related application No. 61/623,876, entitled “IdentificationOf DNA Fragments And Structural Variations” and Ser. No. 13/447,087,entitled “Processing and Analysis of Complex Nucleic Acid SequenceData.”

Definitions

As used herein, the term “complex nucleic acid” refers to largepopulations of nonidentical nucleic acids or polynucleotides. In certainembodiments, the target nucleic acid is genomic DNA; exome DNA (a subsetof whole genomic DNA enriched for transcribed sequences which containsthe set of exons in a genome); a transcriptome (i.e., the set of allmRNA transcripts produced in a cell or population of cells, or cDNAproduced from such mRNA); a methylome (i.e., the population ofmethylated sites and the pattern of methylation in a genome); an exome(i.e., protein-coding regions of a genome selected by an exon capture orenrichment method; a microbiome; a mixture of genomes of differentorganisms; a mixture of genomes of different cell types of an organism;and other complex nucleic acid mixtures comprising large numbers ofdifferent nucleic acid molecules (examples include, without limitation,a microbiome, a xenograft, a solid tumor biopsy comprising both normaland tumor cells, etc.), including subsets of the aforementioned types ofcomplex nucleic acids. In one embodiment, such a complex nucleic acidhas a complete sequence comprising at least one gigabase (Gb) (a diploidhuman genome comprises approximately 6 Gb of sequence).

Nonlimiting examples of complex nucleic acids include “circulatingnucleic acids” (CNA), which are nucleic acids circulating in human bloodor other body fluids, including but not limited to lymphatic fluid,liquor, ascites, milk, urine, stool and bronchial lavage, for example,and can be distinguished as either cell-free (CF) or cell-associatednucleic acids (reviewed in Pinzani et al., Methods 50:302-307, 2010),e.g., circulating fetal cells in the bloodstream of a expecting mother(see, e.g., Kavanagh et al., J. Chromatol. B 878:1905-1911, 2010) orcirculating tumor cells (CTC) from the bloodstream of a cancer patient(see, e.g., Allard et al., Clin Cancer Res. 10:6897-6904, 2004). Anotherexample is genomic DNA from a single cell or a small number of cells,such as, for example, from biopsies (e.g., fetal cells biopsied from thetrophectoderm of a blastocyst; cancer cells from needle aspiration of asolid tumor; etc.). Another example is pathogens, e.g., bacteria cells,virus, or other pathogens, in a tissue, in blood or other body fluids,etc.

As used herein, the term “target nucleic acid” (or polynucleotide) or“nucleic acid of interest” refers to any nucleic acid (orpolynucleotide) suitable for processing and sequencing by the methodsdescribed herein. The nucleic acid may be single stranded or doublestranded and may include DNA, RNA, or other known nucleic acids. Thetarget nucleic acids may be those of any organism, including but notlimited to viruses, bacteria, yeast, plants, fish, reptiles, amphibians,birds, and mammals (including, without limitation, mice, rats, dogs,cats, goats, sheep, cattle, horses, pigs, rabbits, monkeys and othernon-human primates, and humans). A target nucleic acid may be obtainedfrom an individual or from a multiple individuals (i.e., a population).A sample from which the nucleic acid is obtained may contain a nucleicacids from a mixture of cells or even organisms, such as: a human salivasample that includes human cells and bacterial cells; a mouse xenograftthat includes mouse cells and cells from a transplanted human tumor;etc.

Target nucleic acids may be unamplified or they may be amplified by anysuitable nucleic acid amplification method known in the art. Targetnucleic acids may be purified according to methods known in the art toremove cellular and subcellular contaminants (lipids, proteins,carbohydrates, nucleic acids other than those to be sequenced, etc.), orthey may be unpurified, i.e., include at least some cellular andsubcellular contaminants, including without limitation intact cells thatare disrupted to release their nucleic acids for processing andsequencing. Target nucleic acids can be obtained from any suitablesample using methods known in the art. Such samples include but are notlimited to: tissues, isolated cells or cell cultures, bodily fluids(including, but not limited to, blood, urine, serum, lymph, saliva, analand vaginal secretions, perspiration and semen); air, agricultural,water and soil samples, etc.

High coverage in shotgun sequencing is desired because it can overcomeerrors in base calling and assembly. As used herein, for any givenposition in an assembled sequence, the term “sequence coverageredundancy,” “sequence coverage” or simply “coverage” means the numberof reads representing that position. It can be calculated from thelength of the original genome (G), the number of reads (N), and theaverage read length (L) as N×L/G. Coverage also can be calculateddirectly by making a tally of the bases for each reference position. Fora whole-genome sequence, coverage is expressed as an average for allbases in the assembled sequence. Sequence coverage is the average numberof times a base is read (as described above). It is often expressed as“fold coverage,” for example, as in “40-fold (or 40×) coverage,” meaningthat each base in the final assembled sequence is represented on anaverage of 40 reads.

As used herein, term “call rate” means a comparison of the percent ofbases of the complex nucleic acid that are fully called, commonly withreference to a suitable reference sequence such as, for example, areference genome. Thus, for a whole human genome, the “genome call rate”(or simply “call rate”) is the percent of the bases of the human genomethat are fully called with reference to a whole human genome reference.An “exome call rate” is the percent of the bases of the exome that arefully called with reference to an exome reference. An exome sequence maybe obtained by sequencing portions of a genome that have been enrichedby various known methods that selectively capture genomic regions ofinterest from a DNA sample prior to sequencing. Alternatively, an exomesequence may be obtained by sequencing a whole human genome, whichincludes exome sequences. Thus, a whole human genome sequence may haveboth a “genome call rate” and an “exome call rate.” There is also a “rawread call rate” that reflects the number of bases that get an A/C/G/Tdesignation as opposed to the total number of attempted bases.(Occasionally, the term “coverage” is used in place of “call rate,” butthe meaning will be apparent from the context).

As used herein, the term “haplotype” means a combination of alleles atadjacent locations (loci) on the chromosome that are transmittedtogether or, alternatively, a set of sequence variants on a singlechromosome of a chromosome pair that are statistically associated. Everyhuman individual has two sets of chromosomes, one paternal and the othermaternal. Usually DNA sequencing results only in genotypic information,the sequence of unordered alleles along a segment of DNA. Inferring thehaplotypes for a genotype separates the alleles in each unordered pairinto two separate sequences, each called a haplotype. Haplotypeinformation is necessary for many different types of genetic analysis,including disease association studies and making inference on populationancestries.

As used herein, the term “phasing” (or resolution) means sortingsequence data into the two sets of parental chromosomes or haplotypes.Haplotype phasing refers to the problem of receiving as input a set ofgenotypes for one individual ora population, i.e., more than oneindividual, and outputting a pair of haplotypes for each individual, onebeing paternal and the other maternal. Phasing can involve resolvingsequence data over a region of a genome, or as little as two sequencevariants in a read or contig, which may be referred to as local phasing,or microphasing. It can also involve phasing of longer contigs,generally including greater than about ten sequence variants, or even awhole genome sequence, which may be referred to as “universal phasing.”Optionally, phasing sequence variants takes place during genomeassembly.

As used herein, the term “transposon” or “transposable element” means aDNA sequence that can change its position within the genome. In aclassic transposition reaction, a transposase catalyzes the randominsertion of excised transposons into DNA targets. During cut-and-pastetransposition, a transposase makes random, staggered double-strandedbreaks in the target DNA and covalently attaches the 3′ end of thetransferred transposon strand to the 5′ end of the target DNA. Thetransposase/transposon complex inserts an arbitrary DNA sequence at thepoint of insertion of the transposon into the target nucleic acid.Transposons that insert randomly into the target nucleic acid sequenceare preferred. Several transposons have been described and use in invitro transposition systems. For example, in the Nextera™ technology(Nature Methods 6, Nov. 2009; Epicentre Biotechnologies, Madison, Wis.)The entire complex is not necessary for insertion; free transposon endsare sufficient for integration. When free transposon ends are used, thetarget DNA is fragmented and the transferred strand of the transposonend oligonucleotide is covalently attached to the 5′ end of the targetfragment. The transposon ends can be modified by addition of desiredsequences, such as PCR primer binding sites, bar codes/tags, etc. Thesize distribution of the fragments can be controlled by changing theamounts of transposase and transposon ends. Exploiting transposon endswith appended sequences results in DNA libraries that can be used inhigh-throughput sequencing.

Use of Microdroplets and Emulsions

In some embodiments, the methods of the present invention are performedin emulsion or microfluidic devices.

A reduction of volumes down to picoliter levels can achieve an evengreater reduction in reagent and computational costs. In someembodiments, this level of cost reduction is accomplished through thecombination of the MT process with emulsion or microfluidic-typedevices. The ability to perform all enzymatic steps in the same reactionwithout DNA purification facilitates the ability to miniaturize andautomate this process and results in adaptability to a wide variety ofplatforms and sample preparation methods.

Recent studies have also suggested an improvement in GC bias afteramplification (e.g., by MDA) and a reduction in background amplificationby decreasing the reaction volumes down to nanoliter size.

There are currently several types of microfluidics devices (e.g.,devices sold by Advanced Liquid Logic, Morrisville, N.C.) orpico/nano-droplet (e.g., RainDance Technologies, Lexington, Mass.) thathave pico-/nano-drop making, fusing (3000/second) and collectingfunctions and could be used in such embodiments of MT.

Amplifying

According to one embodiment, the MT process begins with a shorttreatment of genomic DNA with a 5′ exonuclease to create 3′single-stranded overhangs that serve as MDA initiation sites. The use ofthe exonuclease eliminates the need for a heat or alkaline denaturationstep prior to amplification without introducing bias into the populationof fragments. Alkaline denaturation can be combined with the 5′exonuclease treatment, which results in a further reduction in bias. Thefragments are amplified, e.g., using an MDA method. In certainembodiments, the MDA reaction is a modified phi29 polymerase-basedamplification reaction, although another known amplification method canbe used.

In some embodiments, the MDA reaction is designed to introduce uracilsinto the amplification products. In some embodiments, a standard MDAreaction utilizing random hexamers is used to amplify the fragments ineach well. In many embodiments, rather than the random hexamers, random8-mer primers are used to reduce amplification bias in the population offragments. In further embodiments, several different enzymes can also beadded to the MDA reaction to reduce the bias of the amplification. Forexample, low concentrations of non-processive 5′ exonucleases and/orsingle-stranded binding proteins can be used to create binding sites forthe 8-mers. Chemical agents such as betaine, DMSO, and trehalose canalso be used to reduce bias through similar mechanisms.

Fragmentation

According to one embodiment, after DNA amplification of DNA, theamplification product, or amplicons, is subjected to a round offragmentation. In some embodiments the CoRE method is used to furtherfragment the fragments in each well following amplification. In order touse the CoRE method, the MDA reaction used to amplify the fragments ineach well is designed to incorporate uracils into the MDA products. Thefragmenting of the MDA products can also be achieved via sonication orenzymatic treatment.

If a CoRE method is used to fragment the MDA products, amplified DNA istreated with a mix of uracil DNA glycosylase (UDG), DNAglycosylase-lyase endonuclease VIII, and T4 polynucleotide kinase toexcise the uracil bases and create single base gaps with functional 5′phosphate and 3′ hydroxyl groups. Nick translation through use of apolymerase such as Taq polymerase results in double-stranded blunt endbreaks, resulting in ligatable fragments of a size range dependent onthe concentration of dUTP added in the MDA reaction. In someembodiments, the CoRE method used involves removing uracils bypolymerization and strand displacement by phi29.

Following fragmentation of the MDA products, the ends of the resultantfragments can be repaired. Such repairs can be necessary, because manyfragmentation techniques can result in termini with overhanging ends andtermini with functional groups that are not useful in later ligationreactions, such as 3′ and 5′ hydroxyl groups and/or 3′ and 5′ phosphategroups. In many aspects of the present invention, it is useful to havefragments that are repaired to have blunt ends, and in some cases, itcan be desirable to alter the chemistry of the termini such that thecorrect orientation of phosphate and hydroxyl groups is not present,thus preventing “polymerization” of the target sequences. The controlover the chemistry of the termini can be provided using methods known inthe art. For example, in some circumstances, the use of phosphataseeliminates all the phosphate groups, such that all ends contain hydroxylgroups. Each end can then be selectively altered to allow ligationbetween the desired components. One end of the fragments can then be“activated”, in some embodiments by treatment with alkaline phosphatase.

Mt Using One of a Small Number of Cells as the Source of Complex NucleicAcids

According to one embodiment, an MT method is used to analyze the genomeof an individual cell or a small number of cells (or a similar number ofnuclei isolated from cells). The process for isolating DNA in this caseis similar to the methods described above, but may occur in a smallervolume.

As discussed above, isolating long fragments of genomic nucleic acidfrom a cell can be accomplished by a number of different methods. In oneembodiment, cells are lysed and the intact nucleic are pelleted with agentle centrifugation step. The genomic DNA is then released throughproteinase K and RNase digestion for several hours. The material canthen in some embodiments be treated to lower the concentration ofremaining cellular waste—such treatments are well known in the art andcan include without limitation dialysis for a period of time (e.g., from2-16 hours) and/or dilution. Since such methods of isolating the nucleicacid does not involve many disruptive processes (such as ethanolprecipitation, centrifugation, and vortexing), the genomic nucleic acidremains largely intact, yielding a majority of fragments that havelengths in excess of 150 kilobases. In some embodiments, the fragmentsare from about 100 to about 750 kilobases in lengths. In furtherembodiments, the fragments are from about 150 to about 600, about 200 toabout 500, about 250 to about 400, and about 300 to about 350 kilobasesin length.

Once isolated, the genomic DNA must be carefully fragmented to avoidloss of material, particularly to avoid loss of sequence from the endsof each fragment, since loss of such material will result in gaps in thefinal genome assembly. In some cases, sequence loss is avoided throughuse of an infrequent nicking enzyme, which creates starting sites for apolymerase, such as phi29 polymerase, at distances of approximately 100kb from each other. As the polymerase creates the new DNA strand, itdisplaces the old strand, with the end result being that there areoverlapping sequences near the sites of polymerase initiation, resultingin very few deletions of sequence.

In some embodiments, a controlled use of a 5′ exonuclease (either beforeor during the MDA reaction) can promote multiple replications of theoriginal DNA from the single cell and thus minimize propagation of earlyerrors through copying of copies.

In one aspect, methods of the present invention produce quality genomicdata from single cells. Assuming no loss of DNA, there is a benefit tostarting with a low number of cells (10 or less) instead of using anequivalent amount of DNA from a large preparation. Starting with lessthan 10 cells ensures uniform coverage in long fragments of any givenregion of the genome. Starting with five or fewer cells allows fourtimes or greater coverage per each 100 kb DNA fragment withoutincreasing the total number of reads above 120 Gb (20 times coverage ofa 6 Gb diploid genome). However, a large number of longer DNA fragments(100 kb or longer) are even more important for sequencing from a fewcells, because for any given sequence there are only as many overlappingfragments as the number of starting cells and the occurrence ofoverlapping fragments from both parental chromosomes can be asubstantial loss of information.

The first step in MT is generally low bias whole genome amplification,which can be of particular use in single cell genomic analysis. Due toDNA strand breaks and DNA losses in handling, even single moleculesequencing methods would likely require some level of DNA amplificationfrom the single cell. The difficulty in sequencing single cells comesfrom attempting to amplify the entire genome. Studies performed onbacteria using MDA have suffered from loss of approximately half of thegenome in the final assembled sequence with a fairly high amount ofvariation in coverage across those sequenced regions. This can partiallybe explained as a result of the initial genomic DNA having nicks andstrand breaks which cannot be replicated at the ends and are thus lostduring the MDA process. MT provides a solution to this problem throughthe creation of long overlapping fragments of the genome prior to MDA.According to one embodiment of the invention, in order to achieve this,a gentle process is used to isolate genomic DNA from the cell. Thelargely intact genomic DNA is then be lightly treated with a frequentnickase, resulting in a semi-randomly nicked genome. Thestrand-displacing ability of phi29 is then used to polymerize from thenicks creating very long (>200 kb) overlapping fragments. Thesefragments are then be used as starting template for MT.

Methylation Analysis Using MT

In a further aspect, methods and compositions of the present inventionare used for genomic methylation analysis. There are several methodscurrently available for global genomic methylation analysis. One methodinvolves bisulfate treatment of genomic DNA and sequencing of repetitiveelements or a fraction of the genome obtained by methylation-specificrestriction enzyme fragmenting. This technique yields information ontotal methylation, but provides no locus-specific data. The next higherlevel of resolution uses DNA arrays and is limited by the number offeatures on the chip. Finally, the highest resolution and the mostexpensive approach requires bisulfate treatment followed by sequencingof the entire genome. Using MT it is possible to sequence all bases ofthe genome and assemble a complete diploid genome with digitalinformation on levels of methylation for every cytosine position in thehuman genome (i.e., 5-base sequencing). Further, MT allow blocks ofmethylated sequence of 100 kb or greater to be linked to sequencehaplotypes, providing methylation haplotyping, information that isimpossible to achieve with any currently available method.

In one non-limiting exemplary embodiment, methylation status is obtainedin a method in which genomic DNA is first denatured for MDA. Next theDNA is treated with bisulfite (a step that requires denatured DNA). Theremaining preparation follows those methods described for example inU.S. application Ser. No. 11/451,692, filed on Jun. 13, 2006 (publishedas US 2007/0072208) and Ser. No. 12/335,168, filed on Dec. 15, 2008(published as US 2009/0311691), each of which is hereby incorporated byreference in its entirety for all purposes and in particular for allteachings related to nucleic acid analysis of mixtures of fragmentsaccording to long fragment read techniques.

In one aspect, MDA will amplify each strand of a specific fragmentindependently yielding for any given cytosine position 50% of the readsas unaffected by bisulfite (i.e., the base opposite of cytosine, aguanine is unaffected by bisulfate) and 50% providing methylationstatus. Reduced DNA complexity helps with accurate mapping and assemblyof the less informative, mostly 3-base (A, T, G) reads.

Bisulfite treatment has been reported to fragment DNA. However, carefultitration of denaturation and bisulfate buffers can avoid excessivefragmenting of genomic DNA. A 50% conversion of cytosine to uracil canbe tolerated in MT allowing a reduction in exposure of the DNA tobisulfite to minimize fragmenting. In some embodiments, some degree offragmenting is acceptable as it would not affect haplotyping.

Using MT for Analysis of Cancer Genomes

It has been suggested that more than 90% of cancers harbor significantlosses or gains in regions of the human genome, termed aneuploidy, withsome individual cancers having been observed to contain in excess offour copies of some chromosomes. This increased complexity in copynumber of chromosomes and regions within chromosomes makes sequencingcancer genomes substantially more difficult. The ability of MTtechniques to sequence and assemble very long (>100 kb) fragments of thegenome makes it well suited for the sequencing of complete cancergenomes.

Error-Reduction by Sequencing a Target Nucleic Acid

According to one embodiment, even if MT-based phasing is not performedand a standard sequencing approach is used, a target nucleic acid isfragmented (if necessary), and the fragments are tagged beforeamplification. An advantage of MT is that errors introduced as a resultof amplification (or other steps) can be identified and corrected bycomparing the sequence obtained from multiple overlapping longfragments. For example, a base call (e.g., identifying a particular basesuch as A, C, G, or T) at a particular position (e.g., with respect to areference) of the sequence data can be accepted as true if the base callis present in sequence data from two or more long fragments (or otherthreshold number), or in a substantial majority of long fragments (e.g.,in at least 51, 60, 70, or 80 percent), where the denominator can berestricted to the fragments having a base call at the particularposition. A base call can include changing one allele of a het orpotential het. A base call at the particular position can be accepted asfalse if it is present in only one long fragment (or other thresholdnumber of long fragments), or in a substantial minority of longfragments (e.g., less than 10, 5, or 3 fragments or as measure with arelative number, such as 20 or 10 percent). The threshold values can bepredetermined or dynamically determined based on the sequencing data. Abase call at the particular position may be converted/accepted as “nocall” if it is not present in a substantial minority and in asubstantial majority of expected fragments (e.g., in 40-60 percent). Insome embodiments and implementations, various parameters may be used(e.g., in distribution, probability, and/or other functions orstatistics) to characterize what may be considered a substantialminority or a substantial majority of fragments. Examples of suchparameters include, without limitation, one or more of: number of basecalls identifying a particular base; coverage or total number of calledbases at a particular position; number and/or identities of distinctfragments that gave rise to sequence data that includes a particularbase call; total number of distinct fragments that gave rise to sequencedata that includes at least one base call at a particular position; thereference base at the particular position; and others. In oneembodiment, a combination of the above parameters for a particular basecall can be input to a function to determine a score (e.g., aprobability) for the particular base call. The scores can be compared toone or more threshold values as part of determining if a base call isaccepted (e.g., above a threshold), in error (e.g., below a threshold),ora no call (e.g., if all of the scores for the base calls are below athreshold). The determination of a base call can be dependent on thescores of the other base calls.

As one basic example, if a base call of A is found in more than 35% (anexample of a score) of the fragments that contain a read for theposition of interest and a base call of C is found in more than 35% ofthese fragments and the other base calls each have a score of less than20%, then the position can be considered a het composed of A and C,possibly subject to other criteria (e.g., a minimum number of fragmentscontaining a read at the position of interest). Thus, each of the scorescan be input into another function (e.g., heuristics, which may usecomparative or fuzzy logic) to provide the final determination of thebase call(s) for the position.

As another example, a specific number of fragments containing a basecall may be used as a threshold. For instance, when analyzing a cancersample, there may be low prevalence somatic mutations. In such a case,the base call may appear in less than 10% of the fragments covering theposition, but the base call may still be considered correct, possiblysubject to other criteria. Thus, various embodiments can use absolutenumbers or relative numbers, or both (e.g., as inputs into comparativeor fuzzy logic). And, such numbers of fragments can be input into afunction (as mentioned above), as well as thresholds corresponding toeach number, and the function can provide a score, which can also becompared to a one or more thresholds to make a final determination as tothe base call at the particular position.

A further example of an error correction function relates to sequencingerrors in raw reads leading to a putative variant call inconsistent withother variant calls and their haplotypes. If 20 reads of variant A arefound in 9 and 8 fragments belonging to respective haplotypes and 7reads of variant G are found in 6 wells (5 or 6 of which are shared withfragments with A-reads), the logic can reject variant G as a sequencingerror because for the diploid genome only one variant can reside at aposition in each haplotype. Variant A is supported with substantiallymore reads, and the G-reads substantially follow fragments of A-readsindicating that they are most likely generate by wrongly reading Ginstead of A. If G reads are almost exclusively in separate fragmentsfrom A, this can indicates that G-reads are wrongly mapped or they comefrom a contaminating DNA.

Identifying Expansions in Regions with Short Tandem Repeats

A short tandem repeat (STR) in DNA is a segment of DNA with a strongperiodic pattern. STRs occur when a pattern of two or more nucleotidesare repeated and the repeated sequences are directly adjacent to eachother; the repeats may be perfect or imperfect, i.e., there may be a fewbase pairs that do not match the periodic motif. The pattern generallyranges in length from 2 to 5 base pairs (bp). STRs typically are locatedin non-coding regions, e.g., in introns. A short tandem repeatpolymorphism (STRP) occurs when homologous STR loci differ in the numberof repeats between individuals. STR analysis is often used fordetermining genetic profiles for forensic purposes. STRs occurring inthe exons of genes may represent hypermutable regions that are linked tohuman disease (Madsen et al, BMC Genomics 9:410, 2008).

In human genomes (and genomes of other organisms) STRs includetrinucleotide repeats, e.g., CTG or CAG repeats. Trinucleotide repeatexpansion, also known as triplet repeat expansion, is caused by slippageduring DNA replication, and is associated with certain diseasescategorized as trinucleotide repeat disorders such as HuntingtonDisease. Generally, the larger the expansion, the more likely it is tocause disease or increase the severity of disease. This property resultsin the characteristic of “anticipation” seen in trinucleotide repeatdisorders, that is, the tendency of age of onset of the disease todecrease and the severity of symptoms to increase through successivegenerations of an affected family due to the expansion of these repeats.Identification of expansions in trinucleotide repeats may be useful foraccurately predict age of onset and disease progression fortrinucleotide repeat disorders.

Expansion of STRs such as trinucleotide repeats can be difficult toidentify using next-generation sequencing methods. Such expansions maynot map and may be missing or underrepresented in libraries. Using MT,it is possible to see a significant drop in sequence coverage in an STRregion. For example, a region with STRs will characteristically have alower level of coverage as compared to regions without such repeats, andthere will be a substantial drop in coverage in that region if there isan expansion of the region, observable in a plot of coverage versusposition in the genome.

For example, if the sequence coverage is about 20 on average, the regionwith the expansion region will have a significant drop, e.g., to 10 ifthe affected haplotype has zero coverage in the expansion region. Thus,a 50% drop would occur. However, if the sequence coverage for the twohaplotypes is compared, the coverage is 10 in the normal haplotype and 0in the affected haplotype, which is a drop of 10 but an overallpercentage drop of 100%. Or, one can analyze the relative amounts, whichis 2:1 (normal vs. coverage in expansion region) for the combinedsequence coverage, but is 10:0 (haplotype 1 vs. haplotype 2), which isinfinity or zero (depending on how the ratio is formed), and thus alarge distinction.

Diagnostic Use of Sequence Data

Sequence data generated using the methods of the present invention areuseful for a wide variety of purposes. According to one embodiment,sequencing methods of the present invention are used to identify asequence variation in a sequence of a complex nucleic acid, e.g., awhole genome sequence, that is informative regarding a characteristic ormedical status of a patient or of an embryo or fetus, such as the sex ofan embryo or fetus or the presence or prognosis of a disease having agenetic component, including, for example, cystic fibrosis, sickle cellanemia, Marfan syndrome, Huntington's disease, and hemochromatosis orvarious cancers, such as breast cancer, for example. According toanother embodiment, the sequencing methods of the present invention areused to provide sequence information beginning with between one and 20cells from a patient (including but not limited to a fetus or an embryo)and assessing a characteristic of the patient on the basis of thesequence.

Cancer Diagnostics

Whole genome sequencing is a valuable tool in assessing the geneticbasis of disease. A number of diseases are known for which there is agenetic basis, e.g., cystic fibrosis,

One application of whole genome sequencing is to understanding cancer.The most significant impact of next-generation sequencing on cancergenomics has been the ability to re-sequence, analyze and compare thematched tumor and normal genomes of a single patient as well as multiplepatient samples of a given cancer type. Using whole genome sequencingthe entire spectrum of sequence variations can be considered, includinggermline susceptibility loci, somatic single nucleotide polymorphisms(SNPs), small insertion and deletion (indel) mutations, copy numbervariations (CNVs) and structural variants (SVs).

In general, the cancer genome is comprised of the patient's germ lineDNA, upon which somatic genomic alterations have been superimposed.Somatic mutations identified by sequencing can be classified either as“driver” or “passenger” mutations. So-called driver mutations are thosethat directly contribute to tumor progression by conferring a growth orsurvival advantage to the cell. Passenger mutations encompass neutralsomatic mutations that have been acquired during errors in celldivision, DNA replication, and repair; these mutations may be acquiredwhile the cell is phenotypically normal, or following evidence of aneoplastic change.

Historically, attempts have been made to elucidate the molecularmechanism of cancer, and several “driver” mutations, or biomarkers, suchas HER2/neu2, have been identified. Based on such genes, therapeuticregimens have been developed to specifically target tumors with knowngenetic alterations. The best defined example of this approach is thetargeting of HER2/neu in breast cancer cells by trastuzumab (Herceptin).Cancers, however, are not simple monogenetic diseases, but are insteadcharacterized by combinations of genetic alterations that can differamong individuals. Consequently, these additional perturbations to thegenome may render some drug regimens ineffective for certainindividuals.

Cancer cells for whole genome sequencing may be obtained from biopsiesof whole tumors (including microbiopsies of a small number of cells),cancer cells isolated from the bloodstream or other body fluids of apatient, or any other source known in the art.

Pre-Implantation Genetic Diagnosis

One application of the methods of the present invention is forpre-implantation genetic diagnosis. About 2 to 3% of babies born havesome type of major birth defect. The risk of some problems, due toabnormal separation of genetic material (chromosomes), increases withthe mother's age. About 50% of the time these types of problems are dueto Down Syndrome, which is a third copy of chromosome 21 (Trisomy 21).The other half result from other types of chromosomal anomalies,including trisomies, point mutations, structural variations, copy numbervariations, etc. Many of these chromosomal problems result in a severelyaffected baby or one which does not survive even to delivery.

In medicine and (clinical) genetics pre-implantation genetic diagnosis(PGD or PIGD) (also known as embryo screening) refers to procedures thatare performed on embryos prior to implantation, sometimes even onoocytes prior to fertilization. PGD can permit parents to avoidselective pregnancy termination. The term pre-implantation geneticscreening (PGS) is used to denote procedures that do not look for aspecific disease but use PGD techniques to identify embryos at risk due,for example, to a genetic condition that could lead to disease.Procedures performed on sex cells before fertilization may instead bereferred to as methods of oocyte selection or sperm selection, althoughthe methods and aims partly overlap with PGD.

Preimplantation genetic profiling (PGP) is a method of assistedreproductive technology to perform selection of embryos that appear tohave the greatest chances for successful pregnancy. When used for womenof advanced maternal age and for patients with repetitive in vitrofertilization (IVF) failure, PGP is mainly carried out as a screeningfor detection of chromosomal abnormalities such as aneuploidy,reciprocal and Robertsonian translocations, and other abnormalities suchas chromosomal inversions or deletions. In addition, PGP can examinegenetic markers for characteristics, including various disease statesThe principle behind the use of PGP is that, since it is known thatnumerical chromosomal abnormalities explain most of the cases ofpregnancy loss, and a large proportion of the human embryos areaneuploid, the selective replacement of euploid embryos should increasethe chances of a successful IVF treatment. Whole-genome sequencingprovides an alternative to such methods of comprehensive chromosomeanalysis methods as array-comparative genomic hybridization (aCGH),quantitative PCR and SNP microarrays. Whole full genome sequencing canprovide information regarding single base changes, insertions,deletions, structural variations and copy number variations, forexample.

As PGD can be performed on cells from different developmental stages,the biopsy procedures vary accordingly. The biopsy can be performed atall preimplantation stages, including but not limited to unfertilizedand fertilized oocytes (for polar bodies, PBs), on day threecleavage-stage embryos (for blastomeres) and on blastocysts (fortrophectoderm cells).

Sequencing Systems and Data Analysis

In some embodiments, sequencing of DNA samples (e.g., such as samplesrepresenting whole human genomes) may be performed by a sequencingsystem. Two examples of sequencing systems are illustrated in FIG. 5.

FIGS. 5A and 5B are block diagrams of example sequencing systems 190that are configured to perform the techniques and/or methods for nucleicacid sequence analysis according to the embodiments described herein. Asequencing system 190 can include or be associated with multiplesubsystems such as, for example, one or more sequencing machines such assequencing machine 191, one or more computer systems such as computersystem 197, and one or more data repositories such as data repository195. In the embodiment illustrated in FIG. 1A, the various subsystems ofsystem 190 may be communicatively connected over one or more networks193, which may include packet-switching or other types of networkinfrastructure devices (e.g., routers, switches, etc.) that areconfigured to facilitate information exchange between remote systems. Inthe embodiment illustrated in FIG. 5B, sequencing system 190 is asequencing device in which the various subsystems (e.g., such assequencing machine(s) 191, computer system(s) 197, and possibly a datarepository 195) are components that are communicatively and/oroperatively coupled and integrated within the sequencing device.

In some operational contexts, data repository 195 and/or computersystem(s) 197 of the embodiments illustrated in FIGS. 5A and 5B may beconfigured within a cloud computing environment 196. In a cloudcomputing environment, the storage devices comprising a data repositoryand/or the computing devices comprising a computer system may beallocated and instantiated for use as a utility and on-demand; thus, thecloud computing environment provides as services the infrastructure(e.g., physical and virtual machines, raw/block storage, firewalls,load-balancers, aggregators, networks, storage clusters, etc.), theplatforms (e.g., a computing device and/or a solution stack that mayinclude an operating system, a programming language executionenvironment, a database server, a web server, an application server,etc.), and the software (e.g., applications, application programminginterfaces or APIs, etc.) necessary to perform any storage-relatedand/or computing tasks.

It is noted that in various embodiments, the techniques described hereincan be performed by various systems and devices that include some or allof the above subsystems and components (e.g., such as sequencingmachines, computer systems, and data repositories) in variousconfigurations and form factors; thus, the example embodiments andconfigurations illustrated in FIGS. 5A and 5B are to be regarded in anillustrative rather than a restrictive sense.

Sequencing machine 191 is configured and operable to receive targetnucleic acids 192 derived from fragments of a biological sample, and toperform sequencing on the target nucleic acids. Any suitable machinethat can perform sequencing may be used, where such machine may usevarious sequencing techniques that include, without limitation,sequencing by hybridization, sequencing by ligation, sequencing bysynthesis, single-molecule sequencing, optical sequence detection,electro-magnetic sequence detection, voltage-change sequence detection,and any other now-known or later-developed technique that is suitablefor generating sequencing reads from DNA. In various embodiments, asequencing machine can sequence the target nucleic acids and cangenerate sequencing reads that may or may not include gaps and that mayor may not be mate-pair (or paired-end) reads. As illustrated in FIGS.5A and 5B, sequencing machine 191 sequences target nucleic acids 192 andobtains sequencing reads 194, which are transmitted for (temporaryand/or persistent) storage to one or more data repositories 195 and/orfor processing by one or more computer systems 197.

Data repository 195 may be implemented on one or more storage devices(e.g., hard disk drives, optical disks, solid-state drives, etc.) thatmay be configured as an array of disks (e.g., such as a SCSI array), astorage cluster, or any other suitable storage device organization. Thestorage device(s) of a data repository can be configured asinternal/integral components of system 190 or as external components(e.g., such as external hard drives or disk arrays) attachable to system190 (e.g., as illustrated in FIG. 5B), and/or may be communicativelyinterconnected in a suitable manner such as, for example, a grid, astorage cluster, a storage area network (SAN), and/or a network attachedstorage (NAS) (e.g., as illustrated in FIG. 5A). In various embodimentsand implementations, a data repository may be implemented on the storagedevices as one or more file systems that store information as files, asone or more databases that store information in data records, and/or asany other suitable data storage organization.

Computer system 197 may include one or more computing devices thatcomprise general purpose processors (e.g., Central Processing Units, orCPUs), memory, and computer logic 199 which, along with configurationdata and/or operating system (OS) software, can perform some or all ofthe techniques and methods described herein, and/or can control theoperation of sequencing machine 191. For example, any of the methodsdescribed herein (e.g., for error correction, haplotype phasing, etc.)can be totally or partially performed by a computing device including aprocessor that can be configured to execute logic 199 for performingvarious steps of the methods. Further, although method steps may bepresented as numbered steps, it is understood that steps of the methodsdescribed herein can be performed at the same time (e.g., in parallel bya cluster of computing devices) or in a different order. Thefunctionalities of computer logic 199 may be implemented as a singleintegrated module (e.g., in an integrated logic) or may be combined intwo or more software modules that may provide some additionalfunctionalities.

In some embodiments, computer system 197 may be a single computingdevice. In other embodiments, computer system 197 may comprise multiplecomputing devices that may be communicatively and/or operativelyinterconnected in a grid, a cluster, or in a cloud computingenvironment. Such multiple computing devices may be configured indifferent form factors such as computing nodes, blades, or any othersuitable hardware configuration. For these reasons, computer system 197in FIGS. 5A and 5B is to be regarded in an illustrative rather than arestrictive sense.

FIG. 6 is a block diagram of an example computing device 200 that can beconfigured to execute instructions for performing variousdata-processing and/or control functionalities as part of sequencingmachine(s) and/or computer system(s).

In FIG. 6, computing device 200 comprises several components that areinterconnected directly or indirectly via one or more system buses suchas bus 275. Such components may include, but are not limited to,keyboard 278, persistent storage device(s) 279 (e.g., such as fixeddisks, solid-state disks, optical disks, and the like), and displayadapter 282 to which one or more display devices (e.g., such as LCDmonitors, flat-panel monitors, plasma screens, and the like) may becoupled. Peripherals and input/output (I/O) devices, which couple to I/Ocontroller 271, can be connected to computing device 200 by any numberof means known in the art including, but not limited to, one or moreserial ports, one or more parallel ports, and one or more universalserial buses (USBs). External interface(s) 281 (which may include anetwork interface card and/or serial ports) can be used to connectcomputing device 200 to a network (e.g., such as the Internet or a localarea network (LAN)). External interface(s) 281 may also include a numberof input interfaces that can receive information from various externaldevices such as, for example, a sequencing machine or any componentthereof. The interconnection via system bus 275 allows one or moreprocessors (e.g., CPUs) 273 to communicate with each connected componentand to execute (and/or control the execution of) instructions fromsystem memory 272 and/or from storage device(s) 279, as well as theexchange of information between various components. System memory 272and/or storage device(s) 279 may be embodied as one or morecomputer-readable non-transitory storage media that store the sequencesof instructions executed by processor(s) 273, as well as other data.Such computer-readable non-transitory storage media include, but is notlimited to, random access memory (RAM), read-only memory (ROM), anelectro-magnetic medium (e.g., such as a hard disk drive, solid-statedrive, thumb drive, floppy disk, etc.), an optical medium such as acompact disk (CD) or digital versatile disk (DVD), flash memory, and thelike. Various data values and other structured or unstructuredinformation can be output from one component or subsystem to anothercomponent or subsystem, can be presented to a user via display adapter282 and a suitable display device, can be sent through externalinterface(s) 281 over a network to a remote device or a remote datarepository, or can be (temporarily and/or permanently) stored on storagedevice(s) 279.

Any of the methods and functionalities performed by computing device 200can be implemented in the form of logic using hardware and/or computersoftware in a modular or integrated manner. As used herein, “logic”refers to a set of instructions which, when executed by one or moreprocessors (e.g., CPUs) of one or more computing devices, are operableto perform one or more functionalities and/or to return data in the formof one or more results or data that is used by other logic elements. Invarious embodiments and implementations, any given logic may beimplemented as one or more software components that are executable byone or more processors (e.g., CPUs), as one or more hardware componentssuch as Application-Specific Integrated Circuits (ASICs) and/orField-Programmable Gate Arrays (FPGAs), or as any combination of one ormore software components and one or more hardware components. Thesoftware component(s) of any particular logic may be implemented,without limitation, as a standalone software application, as a client ina client-server system, as a server in a client-server system, as one ormore software modules, as one or more libraries of functions, and as oneor more static and/or dynamically-linked libraries. During execution,the instructions of any particular logic may be embodied as one or morecomputer processes, threads, fibers, and any other suitable run-timeentities that can be instantiated on the hardware of one or morecomputing devices and can be allocated computing resources that mayinclude, without limitation, memory, CPU time, storage space, andnetwork bandwidth.

Techniques and Algorithms for the MT Process

Basecalling

In some embodiments, data extraction will rely on two types of imagedata: bright-field images to demarcate the positions of all DNBs on asurface, and sets of fluorescence images acquired during each sequencingcycle. Data extraction software can be used to identify all objects withthe bright-field images and then for each such object, the software canbe used to compute an average fluorescence value for each sequencingcycle. For any given cycle, there are four data points, corresponding tothe four images taken at different wavelengths to query whether thatbase is an A, G, C or T. These raw data points (also referred to hereinas “base calls”) are consolidated, yielding a discontinuous sequencingread for each DNB.

A computing device can assemble the population of identified bases toprovide sequence information for the target nucleic acid and/or identifythe presence of particular sequences in the target nucleic acid. Forexample, the computing device may assemble the population of identifiedbases in accordance with the techniques and algorithms described hereinby executing various logic; an example of such logic is software codewritten in any suitable programming language such as Java, C++, Perl,Python, and any other suitable conventional and/or object-orientedprogramming language. When executed in the form of one or more computerprocesses, such logic may read, write, and/or otherwise processstructured and unstructured data that may be stored in variousstructures on persistent storage and/or in volatile memory; examples ofsuch storage structures include, without limitation, files, tables,database records, arrays, lists, vectors, variables, memory and/orprocessor registers, persistent and/or memory data objects instantiatedfrom object-oriented classes, and any other suitable data structures. Insome embodiments, the identified bases are assembled into a completesequence through alignment of overlapping sequences obtained frommultiple sequencing cycles performed on multiple DNBs. As used herein,the term “complete sequence” refers to the sequence of partial or wholegenomes as well as partial or whole target nucleic acids. In furtherembodiments, assembly methods performed by one or more computing devicesor computer logic thereof utilize algorithms that can be used to “piecetogether” overlapping sequences to provide a complete sequence. In stillfurther embodiments, reference tables are used to assist in assemblingthe identified sequences into a complete sequence. A reference table maybe compiled using existing sequencing data on the organism of choice.For example human genome data can be accessed through the NationalCenter for Biotechnology Information at ftp.ncbi.nih.gov/refseq/release,or through the J. Craig Venter Institute at www.jcvi.org/researchhuref/.All or a subset of human genome information can be used to create areference table for particular sequencing queries. In addition, specificreference tables can be constructed from empirical data derived fromspecific populations, including genetic sequence from humans withspecific ethnicities, geographic heritage, religious orculturally-defined populations, as the variation within the human genomemay slant the reference data depending upon the origin of theinformation contained therein.

In any of the embodiments of the invention discussed herein, apopulation of nucleic acid templates and/or DNBs may comprise a numberof target nucleic acids to substantially cover a whole genome or a wholetarget polynucleotide. As used herein, “substantially covers” means thatthe amount of nucleotides (i.e., target sequences) analyzed contains anequivalent of at least two copies of the target polynucleotide, or inanother aspect, at least ten copies, or in another aspect, at leasttwenty copies, or in another aspect, at least 100 copies. Targetpolynucleotides may include DNA fragments, including genomic DNAfragments and cDNA fragments, and RNA fragments. Guidance for the stepof reconstructing target polynucleotide sequences can be found in thefollowing references, which are incorporated by reference: Lander et al,Genomics, 2: 231-239 (1988); Vingron et al, J. Mol. Biol., 235: 1-12(1994); and like references.

In some embodiments, four images, one for each color dye, are generatedfor each queried position of a complex nucleotide that is sequenced. Theposition of each spot in an image and the resulting intensities for eachof the four colors is determined by adjusting for crosstalk between dyesand background intensity. A quantitative model can be fit to theresulting four-dimensional dataset. A base is called for a given spot,with a quality score that reflects how well the four intensities fit themodel.

Basecalling of the four images for each field can be performed inseveral steps by one or more computing devices or computer logicthereof. First, the image intensities are corrected for background usingmodified morphological “image open” operation. Since the locations ofthe DNBs line up with the camera pixel locations, the intensityextraction is done as a simple read-out of pixel intensities from thebackground corrected images. These intensities are then corrected forseveral sources of both optical and biological signal cross-talks, asdescribed below. The corrected intensities are then passed to aprobabilistic model that ultimately produces for each DNB a set of fourprobabilities of the four possible basecall outcomes. Several metricsare then combined to compute the basecall score using pre-fittedlogistic regression.

Intensity correction: Several sources of biological and opticalcross-talks are corrected using linear regression model implemented ascomputer logic that is executed by one or more computing devices. Thelinear regression was preferred over de-convolution methods that arecomputationally more expensive and produced results with similarquality. The sources of optical cross-talks include filter band overlapsbetween the four fluorescent dye spectra, and the lateral cross-talksbetween neighboring DNBs due to light diffraction at their closeproximities. The biological sources of cross-talks include incompletewash of previous cycle, probe synthesis errors and probe “slipping”contaminating signals of neighboring positions, incomplete anchorextension when interrogating “outer” (more distant) bases from anchors.The linear regression is used to determine the part of DNB intensitiesthat can be predicted using intensities of either neighboring DNBs orintensities from previous cycle or other DNB positions. The part of theintensities that can be explained by these sources of cross-talk is thensubtracted from the original extracted intensities. To determine theregression coefficients, the intensities on the left side of the linearregression model need to be composed primarily of only “background”intensities, i.e., intensities of DNBs that would not be called thegiven base for which the regression is being performed. This requirespre-calling step that is done using the original intensities. Once theDNBs that do not have a particular basecall (with reasonable confidence)are selected, a computing device or computer logic thereof performs asimultaneous regression of the cross-talk sources:

I _(background) ^(Base) ≈I _(DNBneighbor1) ^(Base) + . . . +I_(DNBneighborN) ^(Base) +I _(DNB) ^(Base2) +I _(DNB) ^(Base3) +I _(DNB)^(Base4) +I _(DNBpreviousCycle) ^(Base) +I _(DNBotherPosition1)^(Base) + . . . +I _(DNBotherPositionN) ^(Base)+ε

The neighbor DNB cross-talk is corrected both using the aboveregression. Also, each DNB is corrected for its particular neighborhoodusing a linear model involving all neighbors over all available DNBpositions.

Basecall Probabilities:

Calling bases using maximum intensity does not account for the differentshapes of background intensity distributions of the four bases. Toaddress such possible differences, a probabilistic model was developedbased on empirical probability distributions of the backgroundintensities. Once the intensities are corrected, a computing device orcomputer logic thereof pre-calls some DNBs using maximum intensities(DNBs that pass a certain confidence threshold) and uses thesepre-called DNBs to derive the background intensity distributions(distributions of intensities of DNBs that are not called a given base).Upon obtaining such distributions, the computing device can compute foreach DNB a tail probability under that distribution that describes theempirical probability of the intensity being background intensity.Therefore, for each DNB and each of the four intensities, the computingdevice or logic thereof can obtain and store their probabilities ofbeing background (p_(BG) ^(A), p_(BG) ^(C), p_(BG) ^(G), p_(BG) ^(T)).Then the computing device can compute the probabilities of all possiblebasecall outcomes using these probabilities. The possible basecalloutcomes need to describe also spots that can be double or in generalmultiple-occupied or not occupied by a DNB. Combining the computedprobabilities with their prior probabilities (lower prior formultiple-occupied or empty spots) gives rise to the probabilities of the16 possible outcomes:

$p^{A} = {\frac{!{p_{BG}^{A} + {!{p_{BG}^{C} + p_{BG}^{G} + p_{BG}^{T}}}}}{\sum p}*p_{SingleBase}^{prior}}$$p^{AC} = {\frac{!{p_{BG}^{A} + {!{p_{BG}^{C} + p_{BG}^{G} + p_{BG}^{T}}}}}{\sum p}*p_{DoubleOccupied}^{prior}}$$p^{ACG} = {\frac{!{p_{BG}^{A} + {!{p_{BG}^{C} + {!{p_{BG}^{G} + p_{BG}^{T}}}}}}}{\sum p}*p_{TripleOccupied}^{prior}}$$p^{ACGT} = {\frac{!{p_{BG}^{A} + {!{p_{BG}^{C} + {!{p_{BG}^{G} + {!p_{BG}^{T}}}}}}}}{\sum p}*p_{QuadrupleOccupied}^{prior}}$$p^{N} = {\frac{p_{BG}^{A} + p_{BG}^{C} + p_{BG}^{G} + p_{BG}^{T}}{\sum p}*p_{EmptySpot}^{prior}}$

These 16 probabilities can then be combined to obtain a reduced set offour probabilities for the four possible basecalls. That is:

p _(4base) ^(A) =p ^(A)+½(p ^(AC) +p ^(AG) +p ^(AT))+⅓(p ^(ACG) +p^(ACT) +p ^(AGT))+¼p ^(ACGT)+¼p ^(N)

Score Computation:

Logistic regression was used to derive the score computation formula. Acomputing device or computer logic thereof fitted the logisticregression to mapping outcomes of the basecalls using several metrics asinputs. The metrics included probability ratio between the called baseand the next highest base, called base intensity, indicator variable ofthe basecall identity, and metrics describing the overall clusteringquality of the field. All metrics were transformed to be collinear withlog-odds-ratio between concordant and discordant calls. The model wasrefined using cross-validation. The logit function with the finallogistic regression coefficients was used to compute the scores inproduction.

Mapping and Assembly

In further embodiments, read data is encoded in a compact binary formatand includes both a called base and quality score. The quality score iscorrelated with base accuracy. Analysis software logic, includingsequence assembly software, can use the score to determine thecontribution of evidence from individual bases with a read.

Reads may be “gapped” due to the DNB structure. Gap sizes vary (usually+/−1 base) due to the variability inherent in enzyme digestion. Due tothe random-access nature of cPAL, reads may occasionally have an unreadbase (“no-call”) in an otherwise high-quality DNB. Read pairs are mated.

Mapping software logic capable of aligning read data to a referencesequence can be used to map data generated by the sequencing methodsdescribed herein. When executed by one or more computing devices, suchmapping logic will generally be tolerant of small variations from areference sequence, such as those caused by individual genomicvariation, read errors, or unread bases. This property often allowsdirect reconstruction of SNPs. To support assembly of larger variations,including large-scale structural changes or regions of dense variation,each arm of a DNB can be mapped separately, with mate pairingconstraints applied after alignment.

As used herein, the term “sequence variant” or simply “variant” includesany variant, including but not limited to a substitution or replacementof one or more bases; an insertion or deletion of one or more bases(also referred to as an “indel”); inversion; conversion; duplication, orcopy number variation (CNV); trinucleotide repeat expansion; structuralvariation (SV; e.g., intrachromosomal or interchromosomal rearrangement,e.g., a translocation); etc. In a diploid genome, a “heterozygosity” or“het” is two different alleles of a particular gene in a gene pair. Thetwo alleles may be different mutants or a wild type allele paired with amutant. The present methods can also be used in the analysis ofnon-diploid organisms, whether such organisms are haploid/monoploid(N=1, where N=haploid number of chromosomes), or polyploid, oraneuploid.

Assembly of sequence reads can in some embodiments utilize softwarelogic that supports DNB read structure (mated, gapped reads withnon-called bases) to generate a diploid genome assembly that can in someembodiments be leveraged off of sequence information generating MTmethods of the present invention for phasing heterozygote sites.

Methods of the present invention can be used to reconstruct novelsegments not present in a reference sequence. Algorithms utilizing acombination of evidential (Bayesian) reasoning and de Bruijingraph-based algorithms may be used in some embodiments. In someembodiments, statistical models empirically calibrated to each datasetcan be used, allowing all read data to be used without pre-filtering ordata trimming. Large scale structural variations (including withoutlimitation deletions, translocations, and the like) and copy numbervariations can also be detected by leveraging mated reads.

Phasing MT Data

FIG. 7 describes the main steps in the phasing of MT data. These stepsare as follows:

(1) Graph Construction Using MT Data:

One or more computing devices or computer logic thereof generates anundirected graph, where the vertices represent the heterozygous SNPs,and the edges represent the connection between those heterozygous SNPs.The edge is composed of the orientation and the strength of theconnection. The one or more computing devices may store such graph instorage structures include, without limitation, files, tables, databaserecords, arrays, lists, vectors, variables, memory and/or processorregisters, persistent and/or memory data objects instantiated fromobject-oriented classes, and any other suitable temporary and/orpersistent data structures.

(2) Graph Construction Using Mate Pair Data:

Step 2 is similar to step 1, where the connections are made based on themate pair data, as opposed to the MT data. For a connection to be made,a DNB must be found with the two heterozygous SNPs of interest in thesame read (same arm or mate arm).

(3) Graph Combination:

A computing device or computer logic thereof represents of each of theabove graphs is via an N×N sparse matrix, where N is the number ofcandidate heterozygous SNPs on that chromosome. Two nodes can only haveone connection in each of the above methods. Where the two methods arecombined, there may be up to two connections for two nodes. Therefore,the computing device or computer logic thereof may use a selectionalgorithm to select one connection as the connection of choice. Thequality of the mate-pair data is significantly inferior to that of theMT data. Therefore, only the MT-derived connections are used.

(4) Graph Trimming:

A series of heuristics were devised and applied, by a computing device,to stored graph data in order to remove some of the erroneousconnections. More precisely, a node must satisfy the condition of atleast two connections in one direction and one connection in the otherdirection; otherwise, it is eliminated.

(5) Graph Optimization:

A computing device or computer logic thereof optimized the graph bygenerating the minimum-spanning tree (MST). The energy function was setto −|strength|. During this process, where possible, the lower strengthedges get eliminated, due to the competition with the stronger paths.Therefore, MST provides a natural selection for the strongest and mostreliable connections.

(6) Contig Building:

Once the minimum-spanning tree is generated and/or stored incomputer-readable medium, a computing device or logic thereof canre-orient all the nodes with taking one node (here, the first node)constant. This first node is the anchor node. For each of the nodes, thecomputing device then finds the path to the anchor node. The orientationof the test node is the aggregate of the orientations of the edges onthe path.

(7) Universal Phasing:

After the above steps, a computing device or logic thereof phases eachof the contigs that are built in the previous step(s). Here, the resultsof this part are referred to as pre-phased, as opposed to phased,indicating that this is not the final phasing. Since the first node waschosen arbitrarily as the anchor node, the phasing of the whole contigis not necessarily in-line with the parental chromosomes. For universalphasing, a few heterozygous SNPs on the contig for which trioinformation is available are used. These trio heterozygous SNPs are thenused to identify the alignment of the contig. At the end of theuniversal phasing step, all the contigs have been labeled properly andtherefore can be considered as a chromosome-wide contig.

Contig Making

In order to make contigs, for each heterozygous SNP-pair, a computingdevice or computer logic thereof tests two hypotheses: the forwardorientation and reverse orientation. A forward orientation means thatthe two heterozygous SNPs are connected the same way they are originallylisted (initially alphabetically). A reverse orientation means that thetwo heterozygous SNPs are connected in reverse order of their originallisting. FIG. 8 depicts the pairwise analysis of nearby heterozygousSNPs involving the assignment of forward and reverse orientations to aheterozygous SNP-pair.

Each orientation will have a numerical support, showing the validity ofthe corresponding hypothesis. This support is a function of the 16 cellsof the connectivity matrix shown in FIG. 9, which shows an example ofthe selection of a hypothesis, and the assignment of a score to it. Tosimplify the function, the 16 variables are reduced to 3: Energy1,Energy2 and Impurity. Energy 1 and Energy2 are two highest value cellscorresponding to each hypothesis. Impurity is the ratio of the sum ofall the other cells (than the two corresponding to the hypothesis) tothe total sum of the cells in the matrix. The selection between the twohypotheses is done based on the sum of the corresponding cells. Thehypothesis with the higher sum is the winning hypothesis. The followingcalculations are only used to assign the strength of that hypothesis. Astrong hypothesis is the one with a high value for Energy1 and Energy2,and a low value for Impurity.

The three metrics Energy1, Energy2 and Impurity are fed into a fuzzyinference system (FIG. 10), in order to reduce their effects into asingle value—score—between (and including) 0 and 1. The fuzzyinterference system (FIS) is implemented as a computer logic that can beexecuted by one or more computing devices.

The connectivity operation is done for each heterozygous SNP pair thatis within a reasonable distance up to the expected contig length (e.g.,20-50 Kb). FIG. 6 shows graph construction, depicting some exemplaryconnectivities and strengths for three nearby heterozygous SNPs.

The rules of the fuzzy inference engine are defined as follows:

(1) If Energy1 is small and Energy2 is small, then Score is very small.

(2) If Energy1 is medium and Energy2 is small, then Score is small.

(3) If Energy1 is medium and Energy2 is medium, then Score is medium.

(4) If Energy1 is large and Energy2 is small, then Score is medium.

(5) If Energy1 is large and Energy2 is medium, then Score is large.

(6) If Energy1 is large and Energy2 is large, then Score is very large.

(7) If Impurity is small, then Score is large.

(8) If Impurity is medium, then Score is small.

(9) If Impurity is large, then Score is very small.

For each variable, the definition of Small, Medium and Large isdifferent, and is governed by its specific membership functions.

After exposing the fuzzy inference system (FIS) to each variable set,the contribution of the input set on the rules is propagated through thefuzzy logic system, and a single (de-fuzzified) number is generated atthe output—score. This score is limited between 0 and 1, with 1 showingthe highest quality

After the application of the FIS to each node pair, a computing deviceor computer logic thereof constructs a complete graph. FIG. 11 shows anexample of such graph. The nodes are colored according to theorientation of the winning hypothesis. The strength of each connectionis derived from the application of the FIS on the heterozygous SNP pairof interest. Once the preliminary graph is constructed (the top plot ofFIG. 11), the computing device or computer logic thereof optimizes thegraph (the bottom plot of FIG. 11) and reduces it to a tree. Thisoptimization process is done by making a Minimum Spanning Tree (MST)from the original graph. The MST guarantees a unique path from each nodeto any other node.

FIG. 11 shows graph optimization. In this application, the first node oneach contig is used as the anchor node, and all the other nodes areoriented to that node. Depending on the orientation, each hit would haveto either flip or not, in order to match the orientation of the anchornode. FIG. 12 shows the contig alignment process for the given example.At the end of this process, a phased contig is made available.

At this point in the process of phasing, the two haplotypes areseparated. Although it is known that one of these haplotypes comes fromthe Mom and one from the Dad, it is not known exactly which one comesfrom which parent. In the next step of phasing, a computing device orcomputer logic thereof attempts to assign the correct parental label(Mom/Dad) to each haplotype. This process is referred to as theUniversal Phasing. In order to do so, one needs to know the associationof at least a few of the heterozygous SNPs (on the contig) to theparents. This information can be obtained by doing a Trio(Mom-Dad-Child) phasing. Using the trio's sequenced genomes, some lociwith known parental associations are identified—more specifically whenat least one parent is homozygous. These associations are then used bythe computing device or computer logic thereof to assign the correctparental label (Mom/Dad) to the whole contigs, that is, to performparent-assisted universal phasing (FIG. 13).

In order to guarantee high accuracy, the following may be performed: (1)when possible (e.g., in the case of NA19240), acquiring the trioinformation from multiple sources, and using a combination of suchsources; (2) requiring the contigs to include at least two knowntrio-phased loci; (3) eliminating the contigs that have a series oftrio-mismatches in a row (indicating a segmental error); and (4)eliminating the contigs that have a single trio-mismatch at the end ofthe trio loci (indicating a potential segmental error).

FIG. 14 shows natural contig separations. Whether parental data are usedor not, contigs often do not continue naturally beyond a certain point.Reasons for contig separation are: (1) more than usual DNA fragmentationor lack of amplification in certain areas, (2) low heterozygous SNPdensity, (3) poly-N sequence on the reference genome, and (4) DNA repeatregions (prone to mis-mapping).

FIG. 15 shows Universal Phasing. One of the major advantages ofUniversal Phasing is the ability to obtain the full chromosomal“contigs.” This is possible because each contig (after UniversalPhasing) carries haplotypes with the correct parental labels. Therefore,all the contigs that carry the label Mom can be put on the samehaplotype; and a similar operation can be done for Dad's contigs.

Another of the major advantages of the MT process is the ability todramatically increase the accuracy of heterozygous SNP calling. FIG. 16shows two examples of error detection resulting from the use of the MTprocess. The first example is shown in FIG. 16 (left), in which theconnectivity matrix does not support any of the expected hypotheses.This is an indication that one of the heterozygous SNPs is not really aheterozygous SNP. In this example, the A/C heterozygous SNP is inreality a homozygous locus (A/A), which was mislabeled as a heterozygouslocus by the assembler. This error can be identified, and eithereliminated or (in this case) corrected. The second example is shown inFIG. 17 (right), in which the connectivity matrix for this case supportsboth hypotheses at the same time. This is a sign that the heterozygousSNP calls are not real.

A “healthy” heterozygous SNP-connection matrix is one that has only twohigh cells (at the expected heterozygous SNP positions, i.e., not on astraight line). All other possibilities point to potential problems, andcan be either eliminated, or used to make alternate basecalls for theloci of interest.

Another advantage of the MT process is the ability to call heterozygousSNPs with weak supports (e.g., where it was hard to map DNBs due to thebias or mismatch rate). Since the MT process requires an extraconstraint on the heterozygous SNPs, one could reduce the threshold thata heterozygous SNP call requires in a non-MT assembler. FIG. 17demonstrates an example of this case in which a confident heterozygousSNP call could be made despite a small number of reads. In FIG. 17(right) under a normal scenario the low number of supporting reads wouldhave prevented any assembler to confidently call the correspondingheterozygous SNPs. However, since the connectivity matrix is “clean,”one could more confidently assign heterozygous SNP calls to these loci.

Annotating SNPs in Splice Sites

Introns in transcribed RNAs need to be spliced out before they becomemRNA. Information for splicing is embedded within the sequence of theseRNAs, and is consensus based. Mutations in splicing site consensussequence are causes to many human diseases (Faustino and Cooper, GenesDev. 17:419-437, 2011). The majority of splice sites conform to a simpleconsensus at fixed positions around an exon. In this regard, a programwas developed to annotate Splice Site mutations. In this program,consensus splice position models (www.life.umd.edu/labs/mount/RNAinfo)was used. A look-up is performed for a pattern: CAGIG in the 5′-endregion of an exon (“|” denotes the beginning of exon), and MAG|GTRAG inthe 3′-end region of the same exon (“|” denotes the ending of exon).Here M={A,C}, R={A,G}. Further, splicing consensus positions areclassified into two types: type I, where consensus to the model is 100%required; and type II, where consensus to the model is preserved in >50%cases. Presumably, a SNP mutation in a type I position will cause thesplicing to miss, whereas a SNP in a type II position will only decreasethe efficiency of the splicing event.

The program logic for annotating splice site mutations comprises twoparts. In part I, a file containing model positions sequences from theinput reference genome is generated. In part 2, the SNPs from asequencing project are compared to these model positions sequences andreport any type I and type II mutations. The program logic isexon-centric instead of intron-centric (for convenience in parsing thegenome). For a given exon, in its 5′-end we look for the consensus“cAGg” (for positions −3, −2, −1, 0. 0 means the start of exon). Capitalletters means type I positions, and lower-case letters means type IIpositions). In the 3′-end of the exon, a look-up is performed for theconsensus “magGTrag” (for position sequence −3, −2, −1, 0, 1, 2, 3, 4).Exons from the genome release that do not confirm to these requirementsare simply ignored (˜5% of all cases). These exons fall into other minorclasses of splice-site consensus and are not investigated by the programlogic. Any SNP from the genome sequenced is compared to the modelsequence at these genomic positions. Any mismatch in type I will bereported. Mismatch in type II positions are reported if the mutationdeparts from the consensus.

The above program logic detects the majority of bad splice-sitemutations. The bad SNPs that are reported are definitely problematic.But there are many other bad SNPs causing splicing problem that are notdetected by this program. For example, there are many introns within thehuman genome that do not confirm to the above-mentioned consensus. Also,mutations in bifurcation points in the middle of the intron may alsocause splice problem. These splice-site mutations are not reported.

Annotation of SNPs Affecting Transcription Factor Binding Sites (TFBS).

JASPAR models are used for finding TFBSs from the released human genomesequences (either build 36 or build 37). JASPAR Core is a collection of130 TFBS positional frequency data for vertebrates, modeled as matrices(Bryne et al., Nucl. Acids Res. 36:D102-D106, 2008; Sandelin et al.,Nucl. Acids Res. 23:D91-D94, 2004). These models are downloaded from theJASPAR website(http://jaspargenereg.net/cgi-bin/jaspar_db.pl?rm=browse&db=core&tax_group=vertebrates).These models are converted into Position Weight Matrices (PWMs) usingthe following formula: wi=log 2 [(fi+p Ni1/2)/(Ni+Ni1/2)/p], where: fiis the observed frequency for the specific base at position I; Ni is thetotal observations at the position; and p the background frequency forthe current nucleotide, which is defaulted to 0.25 (Wasserman andSandelin, Nature Reviews, Genetics 5:P276-287, 2004). A specificprogram, mast (meme.sdsc.edu/meme/mast-intro.html), is used to searchsequence segments within the genome for TFBS-sites. A program was run toextract TFBS-sites in the reference genome. The outline of steps is asfollows: (i) For each gene with mRNA, extract [−5000, 1000] putativeTFBS-containing regions from the genome, with 0 being the mRNA startinglocation. (ii) Run mast-search of all PWM-models for the putativeTFBS-containing sequences. (iii) Select those hits above a giventhreshold. (iv) For regions with multiple or overlapping hits, selectonly 1-hit, the one with the highest mast-search score.

With the TFBS model-hits from the reference genome generated and/orstored in suitable computer-readable medium, a computing device orcomputer logic thereof can identify SNPs which are located within thehit-region. These SNPs will impact on the model, and a change in thehit-score. A second program was written to compute such changes in thehit-score, as the segment containing the SNP is run twice into the PWMmodel, once for the reference, and the second time for the one with theSNP substitution. A SNP causing the segment hit score to drop more than3 is identified as a bad SNP.

Selection of genes with two bad SNPs. Genes with bad SNPs are classifiedinto two categories: (1) those affecting the AA-sequence transcribed;and (2) those affecting the transcription binding site. For AA-sequenceaffecting, the following SNP subcategories are included:

(1) Nonsense or Nonstop Variations.

These mutations either cause a truncated protein or an extended protein.In either situation, the function of the protein product is eithercompletely lost or less efficient.

(2) Splice Site Variations.

These mutations cause either the splice site for an intron to bedestroyed (for those positions required to be 100% of a certainnucleotide by the model) or severely diminished (for those sitesrequired to be >50% for a certain nucleotide by the model. The SNPcauses the splice-site nucleotide to mutate to another nucleotide thatis below 50% of consensus as predicted by the splice-site consensussequence model). These mutations will likely produce proteins which aretruncated, missing exons, or severely diminishing in protein productquantity.

(3) Polyphen2 Annotation of AA Variations.

For SNPs that cause change in amino-acid sequence of a protein, but notits length, Polyphen2 (Adzhubei et al., Nat. Methods 7:248-249, 2010)was used as the main annotation tool. Polyphen2 annotates the SNP with“benign”, “unknown, “possibly damaging”, and “probably damaging”. Both“possibly damaging” and “probably damaging” were identified as bad SNPs.These category assignments by Polyphen2 are based on structuralpredictions of the Polyphen2 software.

For transcription-binding site mutations the 75% of maxScore of themodels was used based on the reference genome as a screening forTFBS-binding sites. Any model-hit in the region that is <=75% ofmaxScore are removed. For those remaining, if a SNP causes the hit-scoreto drop 3 or more, it is considered as a detrimental SNP.

Two classes of genes are reported. Class 1 genes are those that had atleast 2-bad AA-affecting mutations. These mutations can be all on asingle allele (Class 1.1), or spread on 2 distinct alleles (Class 1.2).Class 2 genes are a superset of the Class 1 set. Class 2 genes are genescontain at least 2-bad SNPs, irrespective it is AA-affecting orTFBS-site affecting. But a requirement is that at least 1 SNP isAA-affecting. Class 2 genes are those either in Class 1, or those thathave 1 detrimental AA-mutation and 1 or more detrimental TFBS-affectingvariations. Class 2.1 means that all these detrimental mutations arefrom a single allele, whereas Class 2.2 means that detrimental SNPs arecoming from two distinct alleles.

The foregoing techniques and algorithms are applicable to methods forsequencing complex nucleic acids, optionally in conjunction with MTprocessing prior to sequencing (MT in combination with sequencing may bereferred to as “MT sequencing”), which are described in detail asfollows. Such methods for sequencing complex nucleic acids may beperformed by one or more computing devices that execute computer logic.An example of such logic is software code written in any suitableprogramming language such as Java, C++, Perl, Python, and any othersuitable conventional and/or object-oriented programming language. Whenexecuted in the form of one or more computer processes, such logic mayread, write, and/or otherwise process structured and unstructured datathat may be stored in various structures on persistent storage and/or involatile memory; examples of such storage structures include, withoutlimitation, files, tables, database records, arrays, lists, vectors,variables, memory and/or processor registers, persistent and/or memorydata objects instantiated from object-oriented classes, and any othersuitable data structures.

Improving Accuracy in Long-Read Sequencing

In DNA sequencing using certain long-read technologies (e.g., nanoporesequencing), long (e.g., 10-100 kb) read lengths are available butgenerally have high false negative and false positive rates. The finalaccuracy of sequence from such long-read technologies can besignificantly enhanced using haplotype information (complete or partialphasing) according to the following general process.

First, a computing device or computer logic thereof aligns reads to eachother. A large number of heterozygous calls are expected to exist in theoverlap. For example, if two to five 100 kb fragments overlap by aminimum of 10%, this results in >10 kb overlap, which could roughlytranslate to 10 heterozygous loci. Alternatively, each long read isaligned to a reference genome, by which a multiple alignment of thereads would be implicitly obtained.

Once the multiple read alignments have been achieved, the overlap regioncan be considered. The fact that the overlap could include a largenumber (e.g., N=10) of het loci can be leveraged to considercombinations of hets. This combinatorial modality results in a largespace (4^(N) or 4̂N; if N=10, then 4^(N)=˜1 million) of possibilities forthe haplotypes. Of all of these 4^(N) points in the N-dimensional space,only two points are expected to contain biologically viable information,i.e., those corresponding to the two haplotypes. In other words, thereis a noise suppression ratio of 4^(N)/2 (here 1e6/2 or ˜500,000). Inreality, much of this 4^(N) space is degenerate, particularly since thesequences are already aligned (and therefore look alike), and alsobecause each locus does not usually carry more than two possible bases(if it is a real het). Consequently, a lower bound for this space isactually 2^(N) (if N=10, then 2^(N)=˜1000). Therefore, the noisesuppression ratio could only be 2^(N)/2 (here 1000/2=500), which isstill quite impressive. As the number of the false positives and falsenegatives grow, the size of the space expands from 2^(N) to 4^(N), whichin turn results in a higher noise suppression ratio. In other words, asthe noise grows, it will automatically be more suppressed. Therefore,the output products are expected to retain only a very small (and ratherconstant) amount of noise, almost independently from the input noise.(The tradeoff is the yield loss in the noisier conditions.) Of course,these suppression ratios are altered if (1) the errors are systematic(or other data idiosyncrasies), (2) the algorithms are not optimal, (3)the overlapping sections are shorter, or (4) the coverage redundancy isless. N is any integer greater than one, such as 2, 3, 5, 10, or more.

The following methodology is useful for increasing the accuracy of thelong-read sequencing methods, which could have a large initial errorrate.

First, a computing device or computer logic thereof aligns a few reads,for instance 5 reads. Assuming reads are ˜100 kb, and the shared overlapis 10%, this results in a 10 kb overlap in the 5 reads or more, such as10-20 reads. Also assume there is a het in every 1 Kb. Therefore, therewould be a total of 10 hets in this common region.

Next, the computing device or computer logic thereof fills in in aportion (e.g., just non-zero elements) or the whole matrix of alpha¹⁰possibilities (where alpha is between 2 and 4) for the above 10candidate hets. In one implementation only 2 out of alpha¹⁰ cells ofthis matrix should be high density (e.g., as measured by a threshold,which can be predetermined or dynamic). These are the cells thatcorrespond to the real hets. These two cells can be consideredsubstantially noise-free centers. The rest should contain mostly 0 andoccasionally 1 memberships, especially if the errors are not systematic.If the errors are systematic, there may be a clustering event (e.g., athird cell that has more than just 0 or 1), which makes the task moredifficult. However, even in this case, the cluster membership for thefalse cluster should be significantly weaker (e.g., as measured by anabsolute or relative amount) than that of the two expected clusters. Thetrade-off in this case is that the starting point should include moremultiple sequences aligned, which relates directly to having longerreads or larger coverage redundancy.

The above step assume that the two viable clusters are observed amongthe overlapped reads. For a large number of false positives, this wouldnot be the case. If this is the case, in the alpha-dimensional space,the expected two clusters will be blurred, i.e., instead of being singlepoints with high density, they will be blurred clusters of M pointsaround the cells of interest, where these cells of interest are thenoise-free centers that are at the center of the cluster. This enablesthe clustering methods to capture the locality of the expected points,despite the fact that the exact sequence is not represented in eachread. A cluster event may also occur when the clusters are blurred (i.e.there could be more than two centers), but in a similar manner asdescribed above, a score (e.g., the total counts for the cells of acluster) can be used to distinguish a weaker cluster from the two realclusters, for a diploid organism. The two real clusters can be used tocreate contigs, as described herein, for various regions, and thecontigs can be matched into two groups to form haplotypes for a largeregion of the complex nucleic acid.

Finally, the computing device or computer logic thereof thepopulation-based (known) haplotypes can be used to increase confidenceand/or to provide extra guidance in finding the actual clusters. A wayto enable this method is to provide each observed haplotype a weight,and to provide a smaller but non-zero value to the unobservedhaplotypes. By doing so, one achieves a bias toward the naturalhaplotypes that have been observed in the population of interest.

Converting Long Reads to Virtual MT

The algorithms that are designed for MT (including the phasingalgorithm) can be used for long reads by assigning a random virtual tag(with uniform distribution) to each of the long fragments. The virtualtag has the benefit of enabling a true uniform distribution for eachcode. MT cannot achieve this level of uniformity due to the differencein the pooling of the codes and the difference in the decodingefficiency of the codes. A ratio of 3:1 (and up to 10:1) can be easilyobserved in the representation of any two codes in MT. However, thevirtual MT process results in a true 1:1 ratio between any two codes.

In view of the foregoing description, according to one aspect of theinvention, methods are provided for determining a sequence of a complexnucleic acid (for example, a whole genome) of one or more organisms,that is, an individual organism or a population of organisms. Suchmethods comprise: (a) receiving at one or more computing devices aplurality of reads of the complex nucleic acid; and (b) producing, withthe computing devices, an assembled sequence of the complex nucleic acidfrom the reads, the assembled sequence comprising less than 1.0, 0.8,0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.08, 0.07, 0.06, 0.05 or 0.04 falsesingle nucleotide variant per megabase at a call rate of 70, 75, 80, 85,90 or 95 percent or greater, wherein the methods are performed by one ormore computing devices. In some aspects, a computer-readablenon-transitory storage medium stores one or more sequences ofinstructions that comprise instructions which, when executed by one ormore computing devices, cause the one or more computing devices toperform the steps of such methods.

According to one embodiment, in which such methods involve haplotypephasing, the method further comprises identifying a plurality ofsequence variants in the assembled sequence and phasing the sequencevariants (e.g., 70, 75, 80, 85, 90, 95 percent or more of the sequencevariants) to produce a phased sequence, i.e., a sequence whereinsequence variants are phased. Such phasing information can be used inthe context of error correction. For example, according to oneembodiment, such methods comprise identifying as an error a sequencevariant that is inconsistent with the phasing of at least two (or threeor more) phased sequence variants.

According to another such embodiment, in such methods the step ofreceiving the plurality of reads of the complex nucleic acid comprises acomputing device and/or a computer logic thereof receiving a pluralityof reads from each of a plurality of long fragments of the complexnucleic acid. Information regarding such fragments is useful forcorrecting errors or for calling a base that otherwise would have been a“no call.” According to one such embodiment, such methods comprise acomputing device and/or a computer logic thereof calling a base at aposition of said assembled sequence on the basis of preliminary basecalls for the position from two or more long fragments. For example,methods may comprise calling a base at a position of said assembledsequence on the basis of preliminary base calls from at least two, atleast three at least four or more than four long fragments. In someembodiments, such methods may comprise identifying a base call as trueif it is present at least two, at least three, at least four longfragments or more than four long fragments. In some embodiments, suchmethods may comprise identifying a base call as true if it is present atleast a majority (or at least 60%, at least 75%, or at least 80%) of thefragments for which a preliminary base call is made for that position inthe assembled sequence. According to another such embodiment, suchmethods comprise a computing device and/or a computer logic thereofidentifying a base call as true if it is present three or more times inreads from two or more long fragments.

According to another such embodiment, the long fragment from which thereads originate is determined by identifying a tag (or unique pattern oftags) that is associated with the fragment. Such tags optionallycomprise an error-correction or error-detection code (e.g., aReed-Solomon error correction code). According to one embodiment of theinvention, upon sequencing a fragment and tag, the resulting readcomprises tag sequence data and fragment sequence data.

According to another embodiment, such methods further comprise: acomputing device and/or a computer logic thereof providing a firstphased sequence of a region of the complex nucleic acid in the regioncomprising a short tandem repeat; a computing device and/or a computerlogic thereof comparing reads (e.g., regular or mate-pair reads) of thefirst phased sequence of the region with reads of a second phasedsequence of the region (e.g., using sequence coverage); and a computingdevice and/or a computer logic thereof identifying an expansion of theshort tandem repeat in one of the first phased sequence or the secondphased sequence based on the comparison.

According to another embodiment, the method further comprises acomputing device and/or a computer logic thereof obtaining genotype datafrom at least one parent of the organism and producing an assembledsequence of the complex nucleic acid from the reads and the genotypedata.

According to another embodiment, the method further comprises acomputing device and/or a computer logic thereof performing steps thatcomprise: aligning a plurality of the reads for a first region of thecomplex nucleic acid, thereby creating an overlap between the alignedreads; identifying N candidate hets within the overlap; clustering thespace of 2^(N) to 4^(N) possibilities or a selected subspace thereof,thereby creating a plurality of clusters; identifying two clusters withthe highest density, each identified cluster comprising a substantiallynoise-free center; and repeating the foregoing steps for one or moreadditional regions of the complex nucleic acid.

According to another embodiment, such methods further comprise providingan amount of the complex nucleic acid, and sequencing the complexnucleic acid to produce the reads.

According to another embodiment, in such methods the complex nucleicacid is selected from the group consisting of a genome, an exome, atranscriptome, a methylome, a mixture of genomes of different organisms,and a mixture of genomes of different cell types of an organism.

According to another aspect of the invention, an assembled human genomesequence is provided that is produced by any of the foregoing methods.For example, one or more computer-readable non-transitory storage mediastores an assembled human genome sequence that is produced by any of theforegoing methods. According to another aspect, a computer-readablenon-transitory storage medium stores one or more sequences ofinstructions that comprise instructions which, when executed by one ormore computing devices, cause the one or more computing devices toperform any, some, or all of the foregoing methods.

According to another aspect of the invention, methods are provided fordetermining a whole human genome sequence, such methods comprising: (a)receiving, at one or more computing devices, a plurality of reads of thegenome; and (b) producing, with the one or more computing devices, anassembled sequence of the genome from the reads comprising less than 600false heterozygous single nucleotide variants per gigabase at a genomecall rate of 70% or greater. According to one embodiment, the assembledsequence of the genome has a genome call rate of 70% or more and anexome call rate of 70% or greater. In some aspects, a computer-readablenon-transitory storage medium stores one or more sequences ofinstructions that comprise instructions which, when executed by one ormore computing devices, cause the one or more computing devices toperform any of the methods of the invention described herein.

According to another aspect of the invention, methods are provided fordetermining a whole human genome sequence, such methods comprising: (a)receiving, at one or more computing devices, a plurality of reads fromeach of a plurality of long fragments, each long fragment comprising oneor more fragments of the genome; and (b) producing, with the one or morecomputing devices, a phased, assembled sequence of the genome from thereads that comprises less than 1000 false single nucleotide variants pergigabase at a genome call rate of 70% or greater. In some aspects, acomputer-readable non-transitory storage medium stores one or moresequences of instructions that comprise instructions which, whenexecuted by one or more computing devices, cause the one or morecomputing devices to perform such methods.

While this invention has been disclosed with reference to specificaspects and embodiments, it is apparent that other embodiments andvariations of this invention may be devised by others skilled in the artwithout departing from the true spirit and scope of the invention.

1. A method of sequencing a target nucleic acid comprising: (a)combining in a single reaction vessel (i) a plurality of long fragmentsof the target nucleic acid, and (ii) a population of polynucleotides,wherein each polynucleotide comprises a tag and a majority of thepolynucleotides comprise a different tag; (b) introducing into amajority of the long fragments tag-containing sequences from saidpopulation of polynucleotides to produced tagged long fragments, whereineach of the tagged long fragments comprises a plurality of thetag-containing sequences at a selected average spacing, and eachtag-containing sequence comprises a tag; (c) producing a plurality ofsubfragments from each tagged long fragment, wherein each subfragmentcomprises one or more tags; (d) sequencing the subfragments to produce aplurality of sequence reads; (e) assign a majority of the sequence readto corresponding long fragments; and (f) assembling the sequence readsto produce an assembled sequence of the target
 2. The method of claim 1wherein producing the subfragments comprises performing an amplificationreaction to produce a plurality of amplicons from each long fragment. 3.The method of claim 2 wherein each amplicon comprises a tag from each ofthe adjacent introduced sequences and a region of the long fragmentbetween the adjacent introduced sequences.
 4. The method of claim 1comprising combining the long fragments with an excess of the populationof tag-containing sequences.
 5. The method of claim 1 comprisingcombining the long fragments with the tag-containing solution underconditions that are suitable for introduction of a single tag-containingsequence into a majority of the long fragments.
 6. The method of claim 1comprising combining the long fragments with the tag-containing solutionunder conditions that are suitable for introduction of a differenttag-containing sequences into a majority of the long fragments.
 7. Themethod of claim 1 wherein the population of tag-containing sequencescomprises a population of beads, wherein each bead comprises multiplecopies of a single tag-containing sequence.
 8. The method of claim 1wherein the population of tag-containing sequences comprises apopulation of concatemers, each concatemer comprising multiple copies ofa single tag-containing sequence.
 9. The method of claim 1 wherein thetag-containing sequences comprise transposon ends, the method comprisingcombining the long fragments and the tag-containing sequences underconditions that are suitable for transposition of the tag-containingsequences into each of the long fragments.
 10. The method of claim 1wherein the tag-containing sequences comprise a hairpin sequence. 11.The method of claim 1 wherein the target nucleic acid is a complexnucleic acid.
 12. The method of claim 1 wherein the target nucleic acidis a genome of an organism.
 13. The method of claim 12 comprisingdetermining a haplotype of the genome.